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ABSTRACT 


The spheroidal theory developed by Vinti for determining the orbit of an 
artificial satellite of an oblate planet is presented in algorithmic form, in 
which empirically derived initial conditions are used to obtain the coordinate 
and velocity components of an unretarded satellite at any time. A differential 
orbit improvement method utilizing observational data is described. This 
method produces a mean set of orbital elements by an iterated least- squares 
fitting of the equations of condition. The results of preliminary applications 
of the orbit generator and differential correction to two artificial satellites of 
the earth, through use of a high- speed digital electronic computer, are shown 
in tabular and graphical form. 
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ORBITAL PREDICTION AND DIFFERENTIAL CORRECTION 
USING VINTI'S SPHEROIDAL THEORY FOR ARTIFICIAL SATELLITES 

by 

Harvey Walden 
Goddard Space Flight Center 


INTRODUCTION 

The spheroidal method for satellite orbits provides a procedure for calculating the orbit of 
any satellite of an oblate planet, when all forces except those of the primary's gravitational field 
are neglected. Determining the effect of the oblateness of a planet on the orbit of a satellite suf- 
ficiently near the planet so that the forces of other bodies may be neglected is one of the central 
problems of satellite astronomy. 

Vinti, in a series of research papers (listed in the bibliography at the conclusion of this report), 
has found a gravitational potential for the exterior of an axially symmetric oblate planet which is 
able to produce an "intermediary reference orbit" accounting for more than 99.5 percent of the 
deviation of the earth's potential from spherical symmetry. The Vinti potential is a very accurate 
approximation for the earth's gravitational potential, which both satisfies Laplace's equation and 
leads to separability of the Hamilton- Jacobi equation in oblate spheroidal coordinates, the most 
appropriate system for an oblate planet. Use of this form for the potential reduces the problem of 
satellite motion to the analytic solution of quartic polynomials and avoids the use of perturbation 
theory entirely in deriving an accurate intermediary orbit. The Vinti potential is actually much 
closer to the empirically accepted one for the earth than any previously used as the starting point 
of a calculation. In the case of the earth, the resulting intermediary orbit reproduces the even 
zonal harmonics exactly through the second and approximately through the fourth. The secular 
solution can be obtained to arbitrarily high order in the second harmonic oblateness parameter, 
and, by means of rapidly converging infinite series, the periodic solution can easily be obtained , 
through second order. The solution holds for all angles of inclination (in the case of equatorial or 
near- equatorial orbits, certain simplifications can be made in the equations) and contains no critical 
inclination or long-periodic terms. For such a reference orbit, error can never accumulate be- 
cause of the exactness of the secular terms. 

This method of solution for unretarded satellite orbits has been adapted for computational 
purposes on a high-speed digital electronic computer primarily by means of the FORTRAN pro- 
gramming language. This paper provides a computational procedure for determining and correct- 
ing an orbit in algorithmic form, adopting algebraic symbols consistent with those in Vinti's 
papers. A summary of preliminary results utilizing observational data from artificial satellites 
is included. 


INPUT PARAMETERS 

The fundamental physical units employed are those of the canonical Vanguard system. In this 
system, the fundamental unit of length is the earth's equatorial radius (taken to be 6378.388 kilo- 
meters), and the fundamental unit of mass is the terrestrial mass (taken to be 5.983x10^"^ kilograms). 
The fundamental unit of time is adjusted so that the Newtonian gravitational constant, G, is set 
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equal to unity; this process yields a value for the Var^ard unit of time of 806.832 seconds. To 
obtain a physical significance for this time, consider a satellite "orbiting"’ the earth at its sur- 
face. This time imit is then seen to be the time required for such a satellite to traverse one 
radian. 

The inertial coordinate system takes the earth’s polar axis as the Z-axis (which is also 
the planetary axis of symmetry and the axis of rotation). The X-Y plane is the equatorial plane, 
with the X-axis pointing toward the vernal equinox (the first point of Aries), the Y-axis orthog- 
onally to the east to form a right-handed system, and the earth’s center of mass at the origin. 

The following constants are required in the computations: 

GM, where G is the Newtonian gravitational constant and M is the earth’s mass. 
From the preceding remarks, it is seen that = 1 in the Vanguard system. 

c = where is the equatorial radius of the earth (unity in the Vanguard sys- 

tem) and J 2 is the coefficient of the second zonal harmonic in the infinite series 
expansion of the earth’s potential. The value of J 2 is approximately 1.0823x10“^ . 

the initial time 

tf, the final time 

At, the time increment used in generating position and velocity components for equal 
time intervals following t. and preceding 

X., Y., z., the initial conditions of position 

X., Y., Z-, the initial conditions of velocity. Note that the set X., y. , z. ,x. , Y. , z. of initial 
conditions is also referred to as the set of injection conditions ^if t.\ the initial 
time, is taken to be the time of injection of the satellite into orbit. 


COORDINATE CONVERSION 

We now compute: 

The square of the magnitude of the position vector: 

r* =X? +Y? +Z2. 

The dot product of the position and velocity vectors: 

r. r. =X,X. +Y.Y. +Z, Z,. 

The oblate spheroidal coordinates (p, v ) and their time derivatives: 
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Then /Oj and tj. are found by extracting square roots, with the condition that the sign of Vi is 
the same as the sign of Z . . 

^ ^ ^ r. (r^ - c^) + 2 c» Z. 

/(rf - c*)2 + 4c* Z* _ ’ 

. 1 r • . (r*-c*) + 2 c*z,z; 

7 ). = — — ~ r. r. + - • - 

2 c% ^ ‘ * + 4 c» Z? J * 

In the above, 0, but if Vi = 0, then Vi = z. 

. . Y, 

$in0. = — . . - , 

r(Pi + (f “ '^i) 

X; 

COS 0 . . 

/(pj + C*) (1 - 7 ^) 

From the above trigonometric relations, we obtain 0. within the limits 0 < < 277 . 

THE JACOBI CONSTANTS OF GENERALIZED MOMENTA 

Compute 

“^1 = ^ + Y? + Z?)-M Pi (pf + c* 77*)"‘ 

a3=XiVi -Yi^i , 

Oj = |(p* + c*i7*)* + a| - 2ttj c* 1?* (1 - 77*)] (1 - 17?)"'^* . 

FACTORING THE QUARTICS: PRIME CONSTANTS 

Compute 



ko = c* Po* . 

A = - 2k„Poy* [l 'ko (3xgy* - 2x* - 8y* + 4)] , 
B = ko p* (1 - yj) [l + k„ y* (4 - xj)] , 
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v-j* ■ 

b2=VF , 

[l "kj, (3xgy2 - 2x2 -8y2 + 4)] ^ 

P = a-‘ [l - Kvl (4 _ x2) - 16 k2y2 (2y2 - 1) 

- •‘o’^yD (*DyD ■ ■ 20yj + 12)] , 


e = ko y| (3x2 - 4) - 16 k2 y2 (2 y2 - 1) 


-k2x2y2 (2x2y| -5x2-28y2 +20) , 

e = [l - 3^ (1 + g)] 

H ( ) = 0, then Vo = ^ and t)" 2 = k^ x 2 (i _ k^ x*). 

If ( a2 - a2 ) ^ 0, then calculate Vo from 


% 


a? - 2 a, c2 
2 _ 2 i_ 


2(a2 - ap 


1 + 


/ 8 a^c 2 (g 2 ^ 

Y (a2-2ajc2)2 


Also, 


MUTUAL CONSTANTS 


If Vo ^ 0> then compute 


a2 - 2 a, c2 
-2 _ ^ ^ 

^ 2(a2 - ap 


1 - /1 + 


8ttj c2 (aj -ap 

2\T 


(aj - 2 0j c2) 


q = Vo^V^ 




A, = (l-e2)X/2p^^^^P 


[(l-e2)>/2] , 


where (x)is the Legendre polynomial with argument x of degree n, and where Rn(^) = 

The infinite series above (and those that follow) is computed by an iterative method, with com- 
putation of terms ceasing when the absolute value of the ratio of successive terms minus unity 
is less than or equal to some preselected tolerance, i.e., computation ceases when 




- 1 


< € 


where e might be 10“^ . Convergence should be attained by consideration of the first several 
terms, in most cases. To increase computational speed, the first term (for n = 2) of the above 
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series may be given explicitly by 


A>? I / ^1 \ 


If 7?o = 0 (corresponding to an orbit in the equatorial plane), compute instead 

A, = (1 - .>)./> p " E.., [(1 - ea,...] ^ 

n-2 

where the first term (for n = 2) is given explicitly by 


If Vo ^ 0, compute 




A, = (1 - p-i y" U] P„ (^) [(1 - , 


\P 

n = 0 


where P„ and R„ are defined as above, and where the first term (for n = 0) of the above series 
is given explicitly by 


(i -e2)i/2 p-i 


If Vo = Oj compute instead 


A = (1 - e2)l/2 p-l Y' M ) R [(1 - e2)l/2] 

Z_J (n! )2 \2p/ " 2 j , 


(n!) 

where the first term (for n = 0) is given explicitly by 


(1 - e2)i/2 p-l 


Then compute 

CD 

A 3 = (1 - e2)i/2 p-3 2^ D„R„,3 [(l-e2)i/2] , 

n=0 

where, for ^ 0, is computed as follows: 



(n an even integer) 



(n an odd integer). 
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If = 0, use instead 


(n an even integer) 


(n an odd integer). 
Then compute 


m=0 


(4m)! 


PJ [(2m)!] 2 W 


(4m + 2)! 


[(2m + 1)!] 2 

V2pj 


B,=i+J.q2+^q4 , 

‘ 2 16 128 


B. = 1 + — + — q4 

2 4^ 64 ^ ’ 


where 


B, = 1 - (1 - T7:2)-1/2 


CD 

E 

m = 2 


- 2m 


yr. =• 


m— 1 

(2m)! y 

22 "> (m!)2 ^ 22" (n!)2 


(2n)' 


All = I (1 - e2)i/2 p-3 e (- 2 bib2 p + bj) 


e2)i/2 p-3 b4 e2 . 


K (bi/bj) < 1, then compute 


A^i = (1 - e2)l/2 p-l e jbi p'l + (3 b] - b2) p“2 
- 1 Wb 2 (l + ^ P'" + I ^2 p-4 (4 + 3e2)] 

^22 = (1 “ e2)‘/2 p-l e2 (3 b] - b^) p"2 


- 1 t>ib2 p"2 + ^ bj P”'' < 
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^23 = I ~ P'* (" W*’! P‘* *^2 P“*) . 


"lie *"2 ®‘* ’ 

Aji = (1 - P’^ e [2 + b, p'l (3 +fej - p'2 (ib| + c*) (4 + 3e^)j , 

A32 = (1 ■ P~® f { 2^2 + (f ®“*)] ’ 

A33 = (1 - p-^ e 3 ^ p-> b, - 1 p -3 (ib3 + c*)j , 

A 34 = -^ (1 - e-* p-3 (|b|+c3) . 

If (b-i/b^ ) ^ 1 (corresponding to a near-equatorial orbit or an equatorial orbit), then compute 
instead: 

A 21 = (1 - e2)i/2 p“ 1 e [bj p"i + (3 ~ b^^) p"4 (i _t?2)2J ^ 

A22 = I (1 - P'* (3 - b|b; 3 ) o'* p--* (1 - 7 , 2)3 ^ 

Aj 3 and A 24 as given above. 

A3, = (1 - e2)i/2 p-3e [2 + (3 + 1 e^) p‘2 (1 - t,*) - (4 + 3 e^) p‘2 , 

A 33 = (1 - e 2 )i /2 p-3 e 2 [j +|c 2 p -2 (1 - rf^y - +|) c* p'^j , 

A33 = I P‘® [I (1 - ^) - ?^] > 

A 34 = -^ (1 - P-* eO c2. 

Now compute 

27JV, = (- 2a,)»/3 (a + b, + A, + c* 7,g Aj B, Bj"!)'* . 
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If (bj/b2 ) < 1, compute 

2771/2 = (^2 " (a + bj -f- Aj + C^ 7 J^ 

If (bj/b^) > 1, compute instead 

2771/2 = [1 -2a^a2^c2(l _ ^2®2^ (a + bj + Aj + c^vl \ Bj 

Then compute 

e' = ae(a + b^)"^ . 

Note that the condition e' < e must be fulfilled. 

THE JACOBI CONSTANTS OF GENERALIZED COORDINATES 

If e / 0, then compute 

Piifij + vl c*) 


sin E- = 


^ ae /(-2aj)(p2 + Apj + B) 
COS = (1 - Pj a“ e"^ . 

From the above, we obtain E. within the limits 0 < E. <2n. 

(1 - sin E^ 


sin V. = 


‘ 1 - e cos E. 


cos v.^ 


cos Ej - e 
1 - e cos E. 


From the above, we obtain v. within the limits 0 < v. < 277. 


If e = 0, then 


and 

If 77^ ^ 0, then compute 


= 0 


E, =0. 


, \ (-°i + 

cos \p. = 7= . - , 

ci7j,/(-2aj)(7j2 - 7 ) 1 ) 


■ I ''i 

Sin i//. = . 

’Jo 


From the above, we obtain i/<. within the limits 0 < >/<j < 2 -n. 
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(1 - 77*)''^* sin 

A 2 — ~ 2 7 

" /I - 7?; sin^ ip. 


COS - 


cos 


/l - 772 sin2 i/;. 


From the above, we obtain y. within the limits 0 < x. < 27 t , 
If t)q = 0, then: i//. = and . 

Now compute 


sin nv. for n = 2 , 3 , 4 , 


sin n0^ for n = 2, 4 . 


If (b^b^ ) < 1, then compute 


/ 3 ^ - (- 2a^)“^'^2 [bjE. + a(E. ^ e s in E^ ) + v. 

+ Aj^ sin V. + Aj 2 sin 2 vJ + c2(a2 - a2)"i/2 ^3 




i (2 + q2) sin 2^/^. + — sin 4 ^. 
8 64 


- t. 


/3j = - a 2 (- 2 aj)-l /2 [A^ + A^j sin v. + A^^ sin 2v. 

+ Aj 3 sin 3v. + Aj^ sin 4v. ] + (a® - Vo “2 


— + 3q®) sin 2 ^ sin 4 '/'i 

32 ‘256 


Ba^i 


If (bj /bj ) > 1, then compute instead 


/ 3 j = (- 2aj)-l''2 [bj E. + a (Ej - e sin E^) + Aj v. 


+ Ajj sin Vj + Aj 2 sin 2Vj] + 77* a'Ml - 2 a^ o.^ (1 - 77^)] 

- i (2 + q^) sin 2>//j + -1- q* sin 40 j - t. , 

8 64 J 
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^2 = “ ^2 [AjV. + A^i sin v. + sin 2 v. 

+ A23 sin 3 v^ + A24 sin 4 v. ] + [1 - 2a^a“2 (1 - 7?^)]“^''^ |B2^i 

— L q2 (4 3 q 2 ) gin 2^//. + ~ 4 r- sin 40. 

32 ^ 256 " 


If 0 and if {h^/h ^ ) < 1, then compute 


^3 = ‘ “ °^3)‘ *^( 


(1 X. 


+ B30i -^2^ sin 20i 


+ (- 2 aj)“^/^ [A 3 V. 


+ A31 sin V. + A32 sin 2 v. + A33 sin 3 v. + A3^ sin 4 v.] . 

If *^0 / 0 and if (b^/b 2 ) > 1, compute instead 

/?3 = 01-^3 a“l [1 - 2a^a-2 c2(l „ 7^2^]“1/2 _ yjlyl/2 (1 _ 7?“2^-1/2 

+ 83 0 j + CL3 (-2aj)"^/2 [A3 v^ + A3J sin V. + A32 sin 2 v. 


+ A33 sin 3 v^ + A3^ sin 4 v^] , 

If 77 ^ =0, compute instead 

- /®2 ( sg " “ 3 ) " “ 3 ^" 2 aj )-*/2 [ A ^ Vj 


+ Aj, sin Vj + Aj 2 sin 2v.] + c^a^(- 2aj)“l/2 [A^ Vj 


+ A3J sin V. + A32 sin 2 v. + A33 sin 3 v. + A34 sin 4 vJ , 


3 

where sgn a = - — - 

1^1 

THE ORBIT GENERATOR OF POSITION ANO VELOCITY COMPONENTS 

In this section, paxameters arise which are tinae-dependent. Initially, the value for time 

t is equal to t., but on subsequent iterations t = t. + n (At), n = 1, 2, 3, Here At is the 

time increment input parameter used in generating position and velocity components for equal 
time intervals following t. and preceding or coincident with the final time, tj . 

li -q^ ^ 0 and if (bj/b^) < 1 , then compute 


= 2 -nv^ (t + / 3 j - c* ^2 a-* 77* Bj B‘*) , 
lAj = 27rvj [t + Aj*(a + b^ + A^)] . 


10 
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Ji ^ ^ and if (b^ /b^ ) > 1, then compute instead 

"Bj(t + /3j) - 77* Bj 


M, = (- 2 aj) 


1/2 


j^a + bj + Aj)B2 + 77^ A2B 


Kjj^ as given above. 

If 77^ = 0 , then compute instead 

M, = (-2aj)i/2 (t +/Si)(a +bj +Aj)-1 , 

'/',=(!- 2aja-*c2)'/2 [y3^ ^ + /5i)(® + bj + Aj)*‘] . 


We now solve the following equation for (m, + Eg). 

M, +Eg - e' sin (M. + Eg) = M, . 

If we let S = M, + Eg , then we can solve this equation (known as Kepler's equation) by use of the 
iterative Newton- Raphson method. 


8 = 8 - (^n- 

" (i - e' COS 8„) 

(8„ - e' sin 8„ - M,)2 (e'sin 8„) 
2(1 - e' cos 8^)^ 

For the initial value, (8 ) „ = m . Iteration ceases when 

' n^n=0 s 



where e is a pre-selected tolerance (e.g., 10 "^ ). Convergence should be attained with several iterations. 
Now use the anomaly connections 

cos v' = (cos 8 - e)(l - e cos 8)"^ , 
sin v' = (1 - (1 - ecos 8)'^ sin 8. 


From the above, we obtain v' within the limits 0 < v' < 277. The angle v' is then placed within 
the same circle of revolution as the angle 8 (which is not taken modulo 2 77). 

Then compute = v' - . 


11 


If Vo 0 and if (t*i Aa ) compute 

'/'o = (- 2 “i)“ (“i - “•1)*''^ % ^ Aj B" 1 Vq , 

Mj = (a +bj)"l [-(Aj+c2t72 AjBjB^'Mvo 

+ ic* (-2aj)‘/2 (0,2 _ a2)-«/2 tj 3 sin (21/I3 +2i/io)] . 

If 7j(, ^ 0 and if (bj /b^) > 1, then compute instead 

r “I 1/2 

i/io = (-2a^)’l/2 [1 _ 2a^a-2 c2 (1 - 7?2)J 

Mj = (a + b^) ^ (A^ + c t]2 

+ lc2 r-2o.j)i/2 l^i ^ 2aja'2 c2 (1 - + 2i//o)| . 

If 7 )^ =0, then compute instead 

^0 = (1 ^ 2aja’2 c2)l/2 A 2 , 

Ml = - (a + b^)-i Aj . 

Continuing, if 77^ / 0 and (b^/b^) < 1, or if 77^ = 0, then compute 

Ej = (1 - e' cos 8)‘ 1 Mj - i e' (1 - e' cos 8)“^ ^2 sin 8 . 

If 7?^ 0 and if (b^/b^) > 1, then compute instead 

Ej = (1 - e' cos 8)" ^ Mj . 

Now use the anomaly connections again 

cos v" = [cos (8 + Ej) - e] [1 - e cos (8 + Ej)] “ ^ , 

s in v" = (1 ~ e2) 1/2 [x - e cos (8 + Ej)] ” ^ sin (8 + E^) . 

From the above, find the angle v" and place it within the same circle of revolution as the 
angle (8 + E^ ). 

Then compute Vj = v" - v' 

If ^0 ^ 0 and if (b^/b^) < 1, then compute 

= (-2aj)-i''2 (a2 -a2)i/2 ^-1 (a^ v, + A^j sin v' 

+ Ajj sin 2v') +-^q* B‘* sin ( 24 >^ + 2i//q) . 
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Ji 7)^ ^ 0 and if {b^/h^) > 1 , then compute instead 

^1 = (“ 2 aj)“i /2 [1 - 2 a" 2 (1 - t? 2)]1/2 g- 1 

+ Aji sin v' + A22 sin 2 v') + -i- q 2 B‘ ^ sin (2 + 2 i/^q) . 

If 7 ?^ = 0, then compute instead 

= (-2a^)'^/2 (1 _ a"2 c^) 1/2 sin v' + A22 sin 2 v' ) . 

Now if (bj/b^) < 1, we continue this procedure one step further to obtain terms of second order, 
as follows. Compute 

M2 = - (a +bj)“^ |aj Vj + Ajj sin v' + Aj2 sin 2 v' 

+ c2 (-2a^)i/2 (0,2 _ a2)~l/2 j ^3 |bj yp^ 

- 1 q2 sin (2 0^ + 2 0^) + -^ q2 sin ( 4 1//^ + 4 ^0) | , 


Let 


E2 = [1 - e' cos (F + Ej)] M2 


E = 8 + Ej + E 2 


and use the anomaly connections once again. 

cos v'" = (cos E - e) (1 - e cos E)“ ^ , 
sin = (1 - e2)i/2 (j cosE)~^ sinE , 

From the above, find the angle v'" and place it within the same circle of revolution as the angle E. 
Then compute = v"* - v" . 

'P 2 = (-2ai)-i/2 (^2 _ ^2)1/2 ^-1 g.j V, +A 2 , V, cos v' 

+ 2 A22 cos 2 v' + A23 sin 3 v' + ^2^ 4 ) 

+ ['/'i COS ( 2 iA^ + 2 i//q) +-|q^ sin ( 2 i/>^ + 2 </.j) q 2 sin + 4 . 


Finally, let 


V =Nl3 + V(, + Vj + Vj , 

'/' = '/’s + *^0 + + '^2 • 
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Now if (bj /b^ ) > 1, we omit computation of ,V 2 , and \p^ . In such case, these terms be- 

come of the third order and hence negligible. 

Instead, we let 

E = 8 + Ej , 

V = M 3 + + Vj , 

'A = 'As + 'Ao + '/'i . 

Continuing, compute 

sin X = (1 - 77J) 1/2 (1 _ 7^2 g^2 1/2 gin «// , 

cosx= (1 ^ t 7 q sin^ cos <// . 

From the above, find the angle y and place it within the same circle of revolution as the 
angle ip . 

If (bj/b^ ) < 1, then compute 

p = (1 + e cos v)” 1 p . 

If (b^b^)^ 1, then compute instead 

p = a (1 -- e cos E) . 

Then, if / 0 and if (bj /b^) < 1, compute 

77 = 7?^ sin i// , 

^ = ^3 + "^3 (^2 - % [(1 “ y 

+ B3 i// + ^ t ?2 sin - c 2 a3 (-2a^)~l/2 (A^ v 

+ Aji sin V + A32 sin 2 v + A33 sin 3 v + A3^ sin 4 v) . 

If T 7 j, ^ 0 and if (bj /b^ ) > 1, compute instead 

V as given above. 

4> = l3^+a^al^ [l - 2 aja -2 c 2 (1 _ r/ 2 )] (1 - -^2 V 

+ 63!/'] - (- 2 aj)*i/^ (AjV + Ajj sin V 

+ A32 sin 2 V + A33 sin 3 v + Aj^ sin 4 v) . 
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I 


If '^0 “ 0 , compute instead 


7 ? = 0 , 


= +^2 (®Sn ag) + ttj (-2aj)-l/2 (A^ v 

+ Ajj sin V + Ajj sin 2v) - (- 2 a ^)~ (Aj v 

+ A3j sin V + A32 sin 2v + A33 sin 3 v + A3^ sin 4 v) , 

a 

where sgn a, = 

1^3! ■ 

The oblate spheroidal coordinates must satisfy the following conditions: 

P > 0 . 


- 1 < 77 < + 1 . 

Now the coordinates and velocities may be found: 

(i - 'n^) ’ 

p = a e (- 2 a^) ^p2 ^ ^ 1/2 ^p2 ^ ^2 c 2 )' 1 sin E , 

Y = c Vq (- 2 a ^) 1/2 (-7^2 _ ^2^ 1/2 ^p2 ^7^2 ^2)“ 1 cos \p , 

X = /{P^\ c2) (1 - 7)^) COS 0 , 

Y = + c^) () - sin 0 , 

Z = p 77 , 

X = X r_P^ , ^1 -y/M , 

Y = Y \_^p_ _ + X , 

[_p 2^^2 i^^ 2 j 

Z = P 77 + 77 p . 

This completes the algorithm for predicting orthogonal position and velocity components of the 
satellite based upon a set of initial conditions. 
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COMPUTATION OF DIRECTION COSINES 


Often the set of initial conditions (position and velocity components) provided is only ap- 
proximate at best, and thus the orbit that is predicted based on these initial conditions will 
similarly contain inaccuracies. In order to remove these inaccuracies and to account for the 
effects of forces not considered by the analytical theory (e.g., aerodynamic drag, solar radia- 
tion, meteoric bombardment, i etc.), the orbital parameters are corrected by comparison with 
those found by direct observation. The orbit improvement method produces a mean set of 
orbital elements through an iterated least- squares fitting of the differential solution to numer- 
ous observational values. 

To perform the differential correction process, the following data must be available in 
addition to the constants listed in the section on input parameters above: 

f s 1 “ Tp /r^ , the flattening coefficient of the earth (where is the polar radius of the 
earth, and is the equatorial radius), taken to be approximately 1/298.3 = 3.3523299 x 10 

, the rotational rate of the earth in radians per mean solar hour (taken to be 0.26251614) 

(x^). , i = 1, 2, . . . , s, the geodetic (or geographic) longitude of the terrestrial tracking 
stations in radians, as measured eastward from Greenwich (a negative sign must be prefixed 
if measured westward from Greenwich). We assume that there are s tracking stations report- 
ing observational data used for comparison purposes. 

i = 1, 2, . . . , s, the geodetic (or geographic) latitude of the stations in radians, 
measured as positive north of the Equator and as negative south of the Equator (- 77/2 < < +^ 7 / 2 ) 

(h). , i = 1, 2, . . . , s, the altitude of the stations in feet, measured as positive above sea 
level and as negative below sea level 

(x^)^ , 1, 2, . . . , the angle in radians, measured eastward from the vernal equinox (the 

first point of Aries) to the Greenwich meridian at midnight Greenwich mean time for each day, d , 
during the period that observations are provided. The apparent sidereal time (the hour angle of 
the first point of Aries) at midnight Greenwich mean time for each day throughout the year is 
tabulated in "The American Ephemeris and Nautical Almanac."* 

a reference time preceding or coinciding with the time of the first observation pro- 
vided, which is used as the zero point in determining t , the relative observation time. It may 
be the time of injection if the tracking data include observations made during the first several 
orbits of the satellite. The purpose of determining a relative observation time, t , is to elimi- 
nate any reference to the calendar. 

We now describe the observation data cards, which are effectively input for the differ- 
ential correction scheme. There are several methods of recording satellite tracking data; 
we present here one of the most common methods, referred to as Minitrack observation data. 
(Another method is discussed in Appendix A.) A set of observation data of this type includes the 
following parameters for each recorded spacecraft observation: 

t' , the date and time of observation. As given, t' is a calendar time. We remove any 
dependence on the calendar by determining t = t' - , where t is the relative observation 

time and is a reference calendar time. Then t becomes a time interval, measured in Van- 
guard units of time from the zero point . As mentioned above, is chosen so that for all 
observations t > 0. It is convenient to have the choice of coincide with that corresponding 
to the initial position and velocity conditions X., Y., Z., X. , Y., z. . When this choice is made, 
is known as an initial or epoch time. 


♦Published annually by Nautical Almanac Office, U. S. Naval Observatory, Washington, D. C. 
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k , the code number for the tracking station reporting the observation. Generally, the 
range of k is the set of integers 1, 2, 3, . . . , s 

Ljj, the observed direction cosine in the X-direction 

Mq, the observed direction cosine in the Y-direction 

w andw^, the weighting factors corresponding to observations andM^^ respectively. 
This information is optional; if not provided, then it is assumed that and are each unity. 

The coordinate system employed for the observation data is centered at the tracking 
station on the earth’s surface, with the X-Y plane tangent to the surface. It is a right-handed, 
orthogonal system with the X-axis extending in an easterly direction along the line of latitude, 
the Y-axis extending in a northerly direction along the line of longitude, and the Z-axis normal 
to the surface and pointing toward the geodetic zenith. 

We first compute the so-called ’’auxiliary functions” S and C (refer to the ’’Explanatory 
Supplement to the Astronomical Ephemeris and the American Ephemeris and Nautical Alma- 
nac”*) from the relations 


C = [cos2 + (1 ^ f)2 sin^ ^ 


S = (1 - f)2 C . 


Here and in the following, we eliminate use of the subscript ”i” referring to an individual one 
of the s tracking stations, and assume that the computations given are performed for each re- 
spective station. The value of H is then converted from its input units of feet to units of the 
earth’s equatorial radius (the conversion factor is 4.77865x10"®) so as to conform to the 
canonical Vanguard system of units used throughout (see INPUT PARAMETERS). Then 
the geocentric latitude is given by 



Now the geocentric distance of the observation point (i.e., tracking station), in units of the 
earth’s equatorial radius, is found: 

p = [(S + H)2 sin2 + (C + H)2 cos2 6^] 1 / 2 . 

The angle S , between the vernal equinox and the observation meridian plane, is computed in 
radian measure by the following expression: 


^ = (^o)d + ^ 


Here, (X.^)^ is as defined above with the value chosen (designated by the subscript "d”) 
corresponding to the midnight immediately preceding observation time. Also, AT is the com- 
puted time, in hours, between observation time and midnight preceding observation time. Thus, 
the second term in the expression for S accounts for the fact that the Greenwich meridian ro- 
tates while the vernal equinox remains fixed in inertial space. 

The inertial geocentric coordinates of the observation point are now converted from a 
spherical to a Cartesian system by means of the followii^ equations: 


•Prepared jointly by the Nautical Almanac Offices of the United Kingdom and the United States of America, London: HM Stationery 
Office, 1961. 
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= p cos 6^ cos S , 
= p cos sin S , 
Z^= p sin . 


The angle , measured in the observation latitude plane between the vernal equinox and the 
tracking station’s X- coordinate axis, is given in radian measure by 


5 77- 

8+2 


The local coordinates of the satellite , Z^) can now be determined from the iner- 

tial position vector (X, Y ,Z ) computed by the orbit generator. Here ’’local” refers to co- 
ordinates measured at the tracking station. The orbit generator will produce the position com- 
ponents (x, Y , Z ) at the observation time. Then the local coordinates are given by the matrix 
relation 


Ym 


10 0 
0 s in cos 6^ 

0 -cos 6^ sin 


cos \jj^ sin 0 

-sin \fj^ cos 0 

0 0 1 



Here, the difference of column matrices on the extreme right represents a translation from 
the earth’s center to the tracking station position; the center matrix on the right represents 
a rotation in the latitude plane about the polar axis through an angle of to bring the inertial 
X-axis into coincidence with the station’s X^^-axis; the left matrix represents a rotation in the 
longitude plane about the X„-axis through an angle of (tt/2 - ) to bring the inertial Z-axis into 

coincidence with the station’s Z„-axis. This matrix equation is equivalent to 




sinxp^ 


or. 


, explicitly stated, 

X„ = (X - X^) cos + (Y - Y^) sin 4;^ 


, =-(X-X^): 




sin sin 0^ + 

-(Y 


(Y - Y^) cos sin 0^ + (Z 
- Y^) cos 0^ cos + 


- Z^) cos 


(Z - Z^) sin 0^ 


We now find the computed values of the direction cosines (in the X-direction) and ( in the 
Y- direction) in terms of the local coordinates. 




Y., 




2n1/2 
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Of course, the computed value of the third direction cosine, N^,(in the Z-direction) is pre- 
determined by and through the relation 




2n 1/2 


THE STANDARD DEVIATION OF FIT 

The differences between the observed and computed values of the direction cosines can 
now be found: 


AM = Mq . 

These differences are sometimes referred to as "residuals”, although this term is also used 
in a different sense in the method of fitting by least squares. We compute these differences 
for each observation in the set of observation data. The number of observations in the set is 
variable, and may be determined by an input parameter, n. 

The average residual is given by 


R=-i 2^ (ALj +AMj), 

i = 1 

where the subscript "i" ranges over individual observations. 

The standard deviation of the residuals from their mean value is found from 



The standard deviation of the residuals 
and is given by 


[(AL. -R)2 + (AMj -R)2] 

(from zero) is called the standard deviation of fit, 


cr 


f 



[(AL.)2 + (AM.)2] . 


As is customary, the larger multiplicative factor (2n _ 6)'^ is used to indicate the excess 
of simultaneous equations of condition over the number of independent coefficients (see FITTING 
BY METHOD OF LEAST SQUARES). 

We may also determine an acceptable range of values for the residuals, bounded by a 
lower limit, , and an upper limit, , based upon the standard deviation. If a residual falls 
outside this range, it may be rejected, with statistical validity, from inclusion in the fittii^ 
process. For example, we may choose 


= R - i o- , 
= R + i cr . 
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For normal (Gaussian) distributions, 68.27 percent of the cases are included within one stand- 
ard deviation on either side of the mean ( j = 1 above), 95.45 percent of the cases are included 
within two standard deviations ( j = 2 above), and 99.73 percent of the cases are included within 
three standard deviations ( j = 3 above). For moderately skewed distributions, the above per- 
centages may hold approximately. K certain of the residuals are rejected on this statistical 
basis, the standard deviation of the accepted residuals only may be computed as a "working" 
standard deviation of fit. Its value is computed exactly as is cTf above, with certain terms 
omitted in the summation, and should be substantially smaller in magnitude than cr^. 


ANALYTICAL PROCEDURE OF DIFFERENTIAL CORRECTION 


The first-order Taylor series expansion of the equations of condition may be written 


AL : 


Lo - k = 


0 

E 


BL 

c 


Aq^. 


I 


AM = Mp - 


E 



i-1 

where ( i = 1, 2, , . . , 6) are the mean or Izsak elements given below. 

Qj = the semi- major axis. 
q 2 = e, the eccentricity. 

Qg = T 7 o = sin I, where I is the inclination of the orbital plane to the Equator. 

= /^i , corresponds to the negative of the time of passage through perigee in Keplerian 
motion. 


^2 i corresponds to the argument of perigee in Keplerian motion. 

= 13^, corresponds to the right ascension of the ascending node in Keplerian motion. 
We may expand the above partial derivatives by the chain rule as follows: 

3qi -dq. Bq^ Bq. 


Bqj BX„ Bqj BY„ BZ„ Bq. 

From the equations for and in terms of the local coordinates given earlier (refer to COMPU- 
TATION OF DIRECTION COSINES), we find directly 

BL 
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BL 


£ V Y rx2 + y2 + 


BL —3/2 

if =- Vm(XS + Y^+Z 5 ) 


bm bl 

ijT ” ” ’ 


3M 


3Y, 


r = (XS + Y5 + zir^'^ - Hi (Xj, + YS + zir^'^ , 


3Zu 


= -Y„Z„(X^+Y2 


Since the coordinates X^,Y^ , z^ and the angles and 6>p are independent of orbital param- 
eters (and merely geodesic functions), we have the matrix relation 






9Y„ 









COS s in 0 

sini//^sin(9jj cos sin cos^j^ 

s in cos - cos cos sin^P 


BX 


Bq, 


by 

Bq. 


BZ 

Bq. 


We find ~ , by substituting 

dq. Bq. dq. 


p - a(l - e cos E) and j] = sin f 


in the relations 


X = /(/>^ + c^) (1 - ? 7 ^) cos 0 > 
Y + c2) (1 - sin0 , 

Z = p7) , 


and then determining BE/^q. , Bi///dq. , and B^/Bq. . Here E and 0 are uniformizing variables, 
analogous respectively to the eccentric anomaly and the argument of latitude in elliptic motion. 
The parameter is the third oblate spheroidal coordinate, the geocentric right ascension. 
The procedure for determining the eighteen partial derivatives BE/Bqj, B^/Bq. , and B0/Bq. is 
a rather lengthy one which is initiated in the next section. 
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Before embarking upon this procedure, the following comments are in order. The ana- 
lytical partial derivatives in the differential correction given in this report correspond to the 
case where ) < 1 only. This excludes equatorial and near- equatorial orbits. More 

specifically, orbits where the inclination, I, is such that 

0 1 I 1 Ic or 180'’ - Ic - 1- 1®^°’ 

where might be as large as 1®54' for an orbit sufficiently close to the earth, are excluded. 

The analytical partial derivatives for the case in which (bj/b^) ^ 1 are, in fact, simpler in 
form, and they are derived in an analogous manner from the equations for this special case 
presented earlier in this report. 

The differential correction process may be carried through to terms of second order or 
it may be simplified to omit terms of purely second order. This is not quite the same as carry- 
ing the process through to terms of first order, since some second-order effects are included 
even when terms of purely second order are dropped. Generally, the speed of computation in 
the differential correction will be increased considerably by neglecting purely second -order 
terms without risking any great loss in precision of the fin^ differential coefficients for the con- 
ditional equations. It should be remembered, however, that even a slight loss in the precision of 
the differential coefficients may necessitate an additional iteration in the least- squares fitting. An 
option to choose the method desired in the differential correction may be provided by inclusion of 
an input variable assuming either one of two values as appropriate for the choice. The terms that 
are to be omitted in the simplified version will be indicated as such in the sections that follow. 


PRIME CONSTANTS II 

The following parameters are utilized extensively throughout the differential correction 
process and must therefore be evaluated beforehand. In many cases, the parameters are those 
that were computed previously by approximation methods, and are here redetermined by more 
accurate expressions: 


p = a (1 - e2) , 

D = (ap “ c^) (ap - c^Vq) + 4 a^ , 


= D + 4 a^ (1 - T}^) , 


A=-2ac2 D"‘ (ap - (1 " vl) , 

B = c2 7?? D' , 



bj = y§ , 

“i = " J M (a + > 


aj = [fi (a + [ap D' D"* - (1 - , 
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- ^o^)r'} (1 - 7,*)*/=* , 
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^ = _d- 2 M 
3a 3a ’ 


3 ir^ _ n- 2 ^ 
3e “■ 3e ’ 


3 D-* ^ _jj-2 

'^■^0 


3b 


= ‘’i \ ^ ^ ® ^ ’ 

da ^ l^a p - ?7 j da J 

^=ac^ (1-7,2) [aI^‘|E Map -0^7,2) , 

3bi frBir^ 1 T 

(®P-c"’7o)-2c2^o D-‘J (1-7,2)-27,oD-i (ap-c2 7,2)| , 


3b, 

37, 


- = i C7,. (D* ir‘)-i/2 Aj-i 9D;.+ p, c(D'D'*)‘/2 

fo 2 \ 3 t,p 3t,o/ 

.'■■(■•a. 


Bai 

^=2 


3a, 1 3b, 

^=-M(a.b,)-_, 


3a, 1 „ 

» =_,i(a + bj)- 2 — , 

37,g 2 3 t,p 


^=lMa.b,)-.]- 


2pD'D-i + ap/rr* ^+D'^5lij 

y 3a 3a y 


[apD'D-‘ - (1 - 7,2)]'*'^^ 


- (a + bj)-> [ap D'D-‘ - (1 _ t,2)]»/2 , 
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— =_ [^.(a + ^aD'D-‘^+ap 




- (a + bj)-‘ [apD'D-1 - (1 - 7]g)] 


^ =i [^(a +bj)-‘]i/2 I ap I^D-* +D' + 2c2t,o [ap D'D-‘ - c* (1 - 


(a + bj)-i [ap D'D-' - ( 1 - tjJ)] ‘/2 , 


- -2 
“2 3a 


2 3e 


_ _ _2 ^ 
3^0 ~ ^■^0 


3a 3a„ 


apD'D-i - c2 (1 - 7,2) 


2[apD'D-l -c2(l -7 j 2)]2 |_ ap D'D* ‘ - c2(1 _ t,2) 


2pD'D-..p D- |5;.D’ 

\ da da 




apD'D-1 - c2 (1 - 7^2) 


2[apD'D’i -c2(l -t? 2)]2 ap -c2 (1 


aD'D-l ^+ap(D'^ ^ + D' ^ 


I 



3^0 


, 2 x 1/2 


f3a. 

r c^T,^ 1 

\3i7o 

apD'D-‘-c2 (1 -77g)_ 


1/2 




2[apD'D'^ -c 2 ( 1 -t; 2)]2 
- 1?0 ^(1 - ^ o )'^ - 



-1/2 - 

^ "apD'D'‘-c*(l -1^) 



BD ' 

”^0 


^■^2 1 -2 -1 


D"M 2D' + a 


^a j 


^^2 1 

Be “ “2 


= - a c-2 77-1 




^^^2 _ 1 _2 -I /n-1 3D' ^ n' ^D-1 

= — a p c“^ 772 ^ 1 D ^ :: + D 


3 i]o 2 


3 % 3 r?o 


Bq _2 ^^2 


3 q _2 


3 i ?2 

Be 


^ '»;■ ( i-q 

3^70 


Bt ?2 

3 '^o 




Bbj 

Ba j ’ 


=(.*b,)-.(.,b.-e””' 


Bb, 


Be / 


Be ^ 

^^0 


Bbi 

- (a + bj)'^ ae 


’ 


Ba 


3 q ,15 3 \|^ 

8 ^ 32 / 3a 


3e \8 ^ 32 ^ j 9e 


3 D '* 
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3,, 41,3^ 


^Vo \8 32 ’ 

^ = ('iq+-Lq3\^ , 

3a \^2 ^ 16 ^ / 3a 


^Jlq._9. qa'i is . 

3e \2^ 16 ^ 3e 


^ J±q.Aq3\ ^ 

^Vo 16^1 


, 3B, 

±_ = - B-2 1 

3a 2 3a 


12 _-- t ,-2 ^ 


3aVbJ-2 


3a Vp 


3bj 

3b, 

- b, b;2 — 2- 

3a 

' 2 3a 

3b, 


- b, br2 — 2- 


1 2 3e 

db^ 


- b, b;2 — 2. 


' " 3^0 

3b, 

- b, p-2 3p 

3a 


3b, 

- b 0-2 

3e 

^2P 3e 


^ ^ = p-i ^ 
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where (x) is the Legendre polynomial with argument x of degree n ^ and P„ (x) is the deriva- 
tive of the Legendre polynomial with respect to the argument. The definition of R is as given 
previously, viz., \ (x) = ( V^). All infinite series are computed by an iterative method, 

with computation of terms ceasing when the absolute value of the ratio of successive terms 
minus unity is less than or equal to some preselected tolerance. 
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+ (1 - e 2 )i /2 p 






[(1 -e^)‘ 


9A, 3p-^ 

^ = Ajp Aje (1 - 


^ (1 - P- (fj 2] " (7) Pn y Kn [(1 - 


f b V^ /b,\ 

\i) "" \i ) "" 


p-^ e J (1 - e^)' 


^-Ij [(1 - P„ P' [(1 - e=)-‘/2] 


Z] " (t) Pn (^j Pn f(l - , 


© 2 . " It/ 


b^V-1 /b. 


— (— 1 
^Vo Vb^y 


f hA " / b, \ 

vtJ "" 


3 A, 


^ a -2 

3 e ~ 2 3 e . 


' ■ * 3^0 


L 



I 1 1 




.1 ^ 

Ba 


„ Bo . / Bb,\ 

.^- 2 p ->( 3 b ;- b 5 , 5 f . 2 p - (^ ib . ^- b , 


. Bo ^‘’i 


fb . P -’ ( l . i.’)( 3 b, 


+ f t> 3 p -3 (4 + 3 e 2 ) 


Be 


= ^1 (1 " ®^)‘‘ + ®‘‘ ■ P‘‘ l^j 

+ (1 - e2)>/3 p-2 e |- bj p-i ^ + ^ ■ 2P'" (3b* - b*) 

^) + f*32P-* li -^2 -^-2bi 

■'^3P-‘ |?)"f‘33p-3e (b3p-l-b,)|, 


+ f b3p-" (4 + 3e3) 


BA, 


^^0 


— - (1 - p~^ e 


9 u -2 

— b« p 
2 ^ ^ 


fBb- / Bb, 3 b,,\ 

(■ * i*’) f> I; * ""k)* 2 ‘‘‘ ' " ' 

{[■P-* ( 3 bJ-b 3 ) 


^^3 _ _ , p-J ^ + i. _-3 g2 /i _ e2sl/2 

3 a ■ 22 P ^ 4 P e (1 e ; 


‘f p-’b,bj-f p-*b;(6..-)]|| 


Q Bb, O 
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1 
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+ 4 - P‘® e 2 (1 -e 2 )i /2 


{ 


-p-‘ (3h\ -hi) +^p-*bj b" 


-|p-3b4(6 + e2)j|j+ (3bj--|p-»b|^ ^ 

+ [- 9p-‘ b, bj +-|p'2 b^ (6 + e2) -bj j 

+ -l-p'3 e (1 - e*)*/* j^3bj - b^ _ 9 p'* bj b| +-|-p"^ b^ (3 + e*)J , 


3 A, 




= — p“^ (1 - 


{( 3 b.-|p-»l) 


[- 


+ 1 - 9 p~^ bj +-~p*2 b| (6 + e^) 


-.] 


'^Vb 


^^3 3^2 

B a B a 


(1 - 77'2)-3/2 7,-3 + 2 2^ "17^77- 


2m- 1 


(where has been given above in the section titled MUTUAL CONSTANTS), 


Be Be 


(1 - -r7~2)-3/2 7,^3 ^ 2 22 n> y^ 77 


- 2m- 1 
2 


3B, 3t7, 


3^0 ^^0 


(1 - 372 ^)"^^^ •>? 2 ® + 2 22 "’■>'m '^2 


- 2m- 1 


/ 3 y 


m = O ' 




( 




3 ^ 7 here y ^u-,] 

!)2 ° / 


^“^0 2^*" (m !)2 ^ (n 

^ = - 3 A, p- U * <1 - p-» ^ R.., [(I - ^ . 


BD 


where BD^/Ba is computed as follows: 
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If n is an even integer, then 




^M¥)E if) 


2 ( i - m ) / b 2 Y”-*p / bj ' 

If; 


FaR 


Ui-nO/bA^- /b,' 


\P/ \bj • 


If n is an odd integer, then 


3a 3a 


= . 2C.-3 |L <-.)‘-<i (^) " p.. 


2 m+ll 


Mv)L * '>(i; 


2(i-m) /b \2"' 


2m+l 1 Ij 


A fh. 

da lb. 


>‘"(i) 


2(i-n) /M'"”''* /M 

\p/ *"’+‘\b 2 y 


cJ A 

— ^ = - A- + e (1 - + (1 - e^) 1/2 p*3 

3 e 3 e 




“1 


p-3 1(1 -e2)-i^[(l-e2)i/2] ‘’'"D„P',2[(l.e2)“i/2]-^(n + 2) [d - e2)i/2] " " ^ ^ ^ [(1 . e^)- 


where has been given above in the section titled MUTUAL CONSTANTS, and where 3D^/3e is 
computed as follows: 

If n is an even integer, then 


^Dn _ 

3 e 3 e 




E 


(i-n,)- 


.( i - n)-l /bA 2 "> /bj 


( b , 


m =0 ' ' \ \ m=o ' 


If n is an odd integer, then 


Be 


BD 




B e 
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(i - m) 




2m + 1 

P 

^ 2 m + 1 



Mt)D- 

m = 0 


iy""’( 2 m + l) 





2m + 1 
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n " 1 


3^0 ’ 


where is computed as follows: 

If n is an even integer, then 


BD„ BD^i 




= 2 p 
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+c 2 AjBj i(a|-a 2 )i/ 2 B;i^ 


/ B B a_ 

,(a^-a2)-x/.A,B- L ^ 


3B;1 r 3b, 3A, 

, (a3 - a|) 1/2 A, -^ - (a^ a|) 1/2 A^ B- 1 + ^ ^ ^ 


, , . 3A, , , / 3B, 3B:1 

+ c2^2b^b;‘^ + c2^2a^ B-i^+B,^ 


(a + bj + Aj + c 2 7)2 Aj Bj B‘i)‘i 


■^''i +‘^^’?oA2BiB;1)-i ](ai-a2)i/2B;i — 


/ Ba Ba \ BB“^ 

+ (a2 - al)- 1/2 A, B; 1 L-^.a^^U (a2 - a|) 1/2 A, ^ 


3b, 3A, 3A- 

-(a2-a2)i/2A,B;i ^ ^ , c2 7,2 B, B' 1 ^ 




+ B-i^+B,-^ (a+bj +A, +c 2 t, 2 a,B,B; 1 )- 


^= 5 ^ 150 * (a +b, +A, +c 2 7)2 A^B, b; 1 )-i |(“|-<^^)^^"b;i^ 


/ 3a, 3a, \ 3B:1 

. (a2 - a 32 )-l /2 A, B- (a, ^ - a3 ^ U (a2 - a|)l/2 A^ ^ 


r3b, 3A, 3A, 

- (“2 - “3) Aj B* 1 ^ + c 2 7)2 Bj B" 1 

2 3 2 2 I 37,^ 3^^ 0 1 2 3^^ 


/ 3B, 3B:1\ 

+ c2 7,2 Aj ^B; 1 — +Bj -^j +2c2 7 ), AjBj B"iJ (a + bj + Aj +c 2 7,2 Aj Bj B^ i)" 1 


(a 2 - a 2 )l /2 7 )-l A2 B;‘ 
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BA 


11 


B a 


= -3A„ p-‘|£-|(i-e2)t/2 p-3 


['"» *"2 If ^ + 2 (bj p - b|) , 


ba. 


” -|p"^ bj (2bjb2 p -b|) ['(l-e*)-!’'^ e2 + 3(l-e*)i/2 p“i e^ - (1 - 

[_ Be 


Be 


■|(1 -e>)i/>p->,b. 




BA 


^ = -|(l-e^)i/3p-3eb. 


, 3b 3b, 

b,p — + 2(b,p-b|)^ 




“”4a-.=>-P-.<b.(^-|p-.b,|£^, 


B a 8 


^.|(l_e3)i/3p-3e2b3 


(1 - e2)i/2 _ e2)'i/2 e2 


BA 


12 3 


Bb^ 


=fa-e^)-P-e3b3^ . 


BA 


23 


= -A23P"'lf +4(1 P'-'e^ b. 


Ba 33*^ 3a 8 


(3bjbjp-i -4b3p-2) 3P 

B a 


Bb. ab ■ 

-^ -^ + 2 (2b| p-‘ -bj) 


3a 


B A<, 


,-i 3p . 1 


3e -A„P-^i^+i(l-e3)3/3p-4e3b^ 


(3bjbjP'» -4b3p-2)|2.bj^ +2(2b3p-> -b,)-p. 


3 b. 


3b„ 


Be 


+ 1 p-4 g2 b3 (b| p“* - bj) 


3 (1 - e3)i/a _ (1 _ e2)->''2 


3A 


23 1 


^ =i(l_ea)3/2p-4e3b^ 


^b, , 3b, 

-b2^.2(2bap-*_b,)^ 


32?oJ 
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^=A 44 .r^- 5 P-.|f), 


BA 


3e 


4 (1 - e 2 )l /2 _ (1 _ e 2 )-i /2 e* 


^. 4 A b- 1 ^. 

'^Vo ^ 3^0 

DIFFERENTIAL CORRECTION: TIME-VARYING PARTIAL DERIVATIVES WITH RESPECT TO 
ENERGY-MOMENTA VARIABLES 

We now compute the following partial derivatives of time- dependent parameters with re- 
spect to the orbital elements a,e , and Vo. We shall later compute the partial derivatives of 
these same parameters with respect to the remaining orbital elements /^i , /S, , and /Sj. 


BM^ ^ 


Ba"^ BB- BB“^ 


TT" V “2 '^0 Bi b;' 


/ ^B, SBrMi 

- -I B- ^ . a- ^ . a- B, , 


BM 

Bt? 




3 a:» BB, 3 B:‘ 

-1 Lx a- In ?_ 


-^1 ^?/32 (BiB-i-^+a;iB-‘^ + a;*B, 


3^0 


-2v, c2/3^a;‘ t7„B,b;1 , 
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[t +/3j +^2 A;;‘ (a +bj +Aj)] 


+ ''2 [-^2 (a +bi + Aj) (al 


1 ^^ 2 * ^ p^-lhl 

Ba ^ Ba 




3 e 


= 277 - 


[‘ +^2^r (a +b, +A,)] 


* "a ^a (a + b, + A,) 




3e 2 3e 




9 '/', 

~^ V , 


^ “ 2 * a;* (a +bj +Aj)] 


r / 9A’* 9a;*\ /9b, 9A,\11 

[^ 2 ( a + b ,. A .) ) J | 




Be’ . p 
— — sin c 
B a 

B a 

\ 3a 

B8 _ 


Be’ . p\ 

Be 

V 9e 

B8 _ 

/3M 

sin 8 

^^0 

U^o 

Cp’ 

< 

© 

1 

^ 2 \ 



sin 8^ (1 - e' cos 8)“*^ 


Ba 


2 ^ 
^ Ba 


^Ms 
B a ’ 


Bv« r “I BM„ 

= (1 _ e2) sin 8 + sin2 8 (sin v*)“^ (1 -e cos8)*^ - — — , 

Be L J 


Bv« PiP 

= (1 - e2) sin 8 (sin v')“i (1 - e cos 8)“^ t- » 
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3a 


= *• -B-.al 

da 2 3a I 

V 2 a 2 a/ ’ 

L\ " 3e ■ ‘"a 3^j (^|-a|)-* - ^a-x ^ . g-i 

^ 3e 2 3 gJ 

* <- 2.,, -/><»... I,... ,.. j,, ^ ^ ^ Sf 


3e " ''0 


, r/ 3^ 3a,\ 

b H a^j <«l - «r ■ - i.;. . B. . „-,] 

M- 2 .,)- 2 = ,., 3 ., A a^' 


3M. 


3 a " "“i <^ + b,)-' 1 + -— -(a + 


" ^%J ’ 


da) ~ + bj)-i ^(Aj + c2 7,2 g.j^ ^ 


3b 


f3A, 


3A, 




* !!i , A B 

3a ^^2 B, — I 


f 0 ^ 7,3 (- 2 a ,)»^2 2 )- V 2 


i . 3a 

^a, sin (2^^ + 2./- ) ' 

da 


+ cos 


(^'P. + 2 ^ 0 ) _ 1 ,2 _ ^27-1 / 3“2 3«3\ ll 

^ ^ 2 ^ b> ( 2 3g -a3 -^j Sin (2^^ + 2^„)J| . 

3 Mj ^ 

-^ = - M, (a + b,)-> ^ _ (3 , g^)-i L ^ 9^ 

3e 

fe.c.,. (b,b,.^,a.b,.^.a.b.!|:)] 

^ 73 ,-7,,.2. <a. . [. .., ,, ^ ^ ^ ^ ^ 

\ 3e 3e/ 

|<7|-.|.- (s^-a.^) aa„(a..,2a, 


+ V„ 
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3Mj 


Bb, 


-I "i -1 

Mj (a + bj) (a + bj) 


(Aj + Vo A: 






'^Vo 


pA, 

b^o 


+ vl 


Bj BJ‘ 


+ A, Br' + A 

3^0 ^^0 


2B, 


3B-^ 


+ 2c^Vo>^ B, 


i BJ' 


,2 „3 


nS(-2aj) 


1/2 


(a*- 


,-1/2 




Ba- 

1 sin (2^, + 2-Z'o) ^ + cos (2>A, 


+ 2s!-o) 






sin (2<p^ + 2</<o) 


3 2 2/ vl/2 

- -c2 7,2 (« 2a,) 


(a^ - a|) 


- 1/2 


sin (2i//^ + 20o) 


3Ei 

Ba 


(1 - e' cos 8)”^ 


^a 


M, (1 - e' cos 8)' 


^e' sin 


s 


Ba 


- cos 6 



- M, (1 - e' cos Py^ sin P 



3e* 

da 


- -|m, e' (1 - e' cos 8) 


1 


- P 

Sin 8 cos t 

da 



1 . n ^8 , 

+ JL M, e' cot 8 — + e' 

2 ^ oa 


BM, 

Ba 


BEi 

Be 


-1 

(1 - e' cos8) 


Mj (1 - e' cos 8 ) 


sin 8 


Be 



- Mj (1 - e' cos 8)“^ sin 8 



Be' 

Be 


^ e' (1 - e' cos 8 ) 


1 


sin 8 


Be 


- cos 8 



+ 



COt8 


b8 

Be 


BM, 

Be 


— i = (1 e cos 8) 


C21W, _ 

M, (1 - e' cos 8) 


sin 8 


B8 

^^0 


- cos 8 



- Mj (1 - e' cos py^ sin 8 




(1 


e' cos 8) ^ 



- cos 


8 



+ 


4 m , 


e' cot 8 


b8 

^Vo 


» * 
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3v, 

- (1 -e2) sin (8 + Ej) (sin v")‘* 


[1 - e cos (g + Ej)]“^ 



Be 


(1 - e^) sin (8 + Ej) 



+ sin^ (g + Ej) 


(sin v")'^ [1 - e cos (g + Ej)]'^ 



Bvj 

drj^ ~ (£+Ej) (sin v")*' [l-e cos (g+Ej)]'^ 



3i//j 

Ba 


" I q"* Bji sin ( 2 ^^ + 2^o) 



3a 


1 


Ba 


+ 


(“2 - ^|) 


2^-l 




-BJ* 



+ (- 2 aj )-»/2 ( a | 







Ba 


+ (Ajj cos v' + 2 Aj 2 cos 2v') 



+ sin 



+ sin 2v' 



Ba 


.lq^B-cos(2^,.2^„) 

+ ^ q B-J sin (2v!-^ + 2v^„) |3 _ | q2 g -2 (2./-^ + 24 >^) ^ , 
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Be 


^^0 


and 

aMj 

Ba 


1 

1 

Ba. 

/ Ba2 


'Pi “ ■! ®2* + 2<pg) 

a7^ 

2 1 

^ + (“2 - “i) 




BB. 

3'i 1 

2 3e 


+ (a2 - 7,^* B’ ‘ 


Bvj BA2 


I'A ' A ''I I ^'^“1 • ' 

+ (Ajj COS V + 2 A 22 COS 2v ) 1 — + —I + sin V — + sin 2v — 


•I Chxp Bi//A 1 

+ - q2 B-i cos (2./-^ + 2>Ao) + -^) + 5 Q B;» sin (2./.^ + 2-/-o) ^ 


1 o 7 

i q2 B-^ sin (2V-, + 2 ^ 0 ) ^ , 


'/'i “ ■§• ®2 ‘ + 2i/<o) 


1 -1 ^“1 / 2 2x-l ( ^“2 


BB, 

Br* — - - --1 




^0 


+ (-2a^)-'"^ (a2-a2)^/S-lBj> 


Bvj BA2 


+ (A,, cos v' + 2 A 22 cos 2v' ) 


"s “'oi , “"21 . , ^^2 


BM Bv, 




3A, 


+ sin v^ — + sin 2v 


3^0 3t7„ 


^ ('bip 1 

+ i q2 Bjl cos (2^^ + 2V/o) ^ ^ + - q B-> sin {2^^ + 2^^) ^ 


3t]o 37?^/ 4 


3q 


3B. 


1 q2 Bj 2 sin (2V', + 2>Ao) 


The followii^ time- dependent partial derivatives with respect to the orbital elements a, e, 
7 q are computed only if the differential correction is carried through terms of second order: 


M2 (a + bj) 


_ (a + b,) <^A, — + V, — + sin v — 


+ (Aj, cos V 


/BM 

' +2A^, cos 2v') + sin2v' 


3Aj2 

Ba 


c2 773 


<- («i - ^ . <- Cl - -I)-"’ (.. ^ 


®i '/'i 
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- i '/'i cos ( 2 V', + 2 </'o) - j sin ( 2 V<, + 2^^) 


+ ^ sin ( 4 i/', + 4 V'o) 




(- 2 a,) 


1/2 




3 V -1 


Ba 


1 ^'^1 

-cos ( 2 i//, + 2 V'o) + V'l sin ( 2 V', + 2 V'o) 



- iq 2 cos ( 2 V«, + 2 V'o) 



- i q sin ( 2 V', + 2 </'o) ' 1 ^ ( 4 V>, + 4 V'o) 



+ 


q sin (40 

32 “ 


+ 40 o ) 


V 

da 


3 M 2 

3 e 


-1 

^2 (a + b ,) + b ,) 


3 v, 

3 e 


+ V , 


Be 


BA. 


+ sm V 


11 

Be 


+ (Ajj cos v' + 2Aj2 cos 2 v') 



+ sin 2v* 


BA 


12 


Be 


“ 7^3 


(- 2 a ,) 


- 1/2 




Ba 


Be 


' + (- 2 a ,)'/" 


( a ^ - 




- i/^, COS (20^ 


+ 2iAo) - sin (20^ 




o ) 


+ c* vl (- 2 a,) 


1/2 




30 , 


3B, 1 

^ - i cos ( 20 ^ + 20 „) 


30 , 

Be 


+ 0 , sin ( 20 ^ + 20 o) 



+ 



- i q 2 cos ( 20 ^ + 20 (,) 



i q sin ( 20 , + 20 (,) + A q 2 cos ( 40 , + 40 „) 



Be / 


+ ^ q sin ( 4 i//^ + 4 i//q) 


)^U 

0^ 3ejj ’ 
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3Mj 


■ + sin V 




sin 2 v' 


2 a, (a| _ ,|)-3/l ( ^ 

\ ^^0 ^ 3%/ 


" I 3^0 *• 3 t 7 o 

+ (A„ cos v' + 2 A , 2 cos 2 v') 1V\ 

3%/ 

- [^0 (-2a,)->/2 («2 _ ^ 2 yU 2 ^ ^ 

■ 3 (~2aj)*/2 ^ _ cl2^i/21 L . 1 

2 3) j |B^ Kjj^ ~ ly^i cos (2l// + 20 'i ^ n2 O ' 

J L 2 1 ^ - J q sin ( 20 ^ + 20 ^) 

^ sin (4>A, , ^ ^3 L 3 ^ ^ ^ 

L ‘ 37?^ + 'f'l 3^^ - 2 

" f20^ + 2^,) ^ |M 

377 J 

3 ~ + q^ cos (4^^ + 4 ^ ) z'!:^ ^ 

\^% ^VoJ 

+ i q sin (40^ + 40 „) ^j| , 


cos ( 24 >^ + 20„) 


+ 0j sin (20 


4 + 2 '/’o) ^ + -L q2 


SEj 

■^ = [1 - e' cos (8 + E,)] 

* da 


3'/', 


^2 [1 - e' cos (8 + E,)] 


e' sin (8 + E 


■>(l!^5)- 


cos (8 + Ej) ^ 
9a 


^^2 r 
de " ~ 


cos (8 + E,)]"' ^ 
‘ 3e 


Mj [1 - e' cos (8 + E 


i)J ^ fe' sin (8 + E ) 3 il 

L a^j -cos(8.E,) ^J, 


3Ej 

■a_ ~ f 1 - e’ cos I'P J. IT \ 1 ~> 


^Vo 


cos (8 + Ej)]~* i 

3t?o 


- M 2 [l - e' cos (8 


+ * fe' sin (8 + E i (38 ,"| 

L ■’ K'k] 
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Ba Ba 3 a Ba 


Be Be Be Be 


2E.=^+ ^ + ^ 
^-^0 ^"^0 ^'^0 ’ 


Ba 


(1 - e^) sinE(sinv) ^(1-ecosE)" 


Be 

Ba 


[<■ - 


e^) sin E ^5. ^ sin^ e| (sin v)“^ (1 - e cos 
Be J 


\ Ba ^ Ba ^ Ba/ ’ 

' \Be Be Be/ 


Bv 

-= (1 - e^) sin E (sin v)'^ (1 - e cos 

^V„ 


Ba Ba Ba Ba Ba ’ 


^ ^ ^ ^ 

3e 3e ^ Be ^ Be Be 


Bv Bv_ Bv, Bv„ 

= 1 0 ^ ^ 2 

^Vo 3^0 ^'^0 ^^0 ^'^0 


<BM 

S 

3^0 

Bv; 





B'/'j 


|>/'2 - -| ®2^ ['^1 + 2.//g) +|- sin (2<A^ + 2i/'o) 

- ^ + 4^o)]} [- f “V It " - “ 3>-‘ (“2 S - “3 ^) - 5 

+ (- 2 a ^)“^/2 ( q _2 ^ a 2 ) l /2 ^-1 g-i 

L 

/bm^ BvA 

- (A 21 V^ sin v' + 4 A 22 v^ sin 2 v' - 3 A ^3 cos 3 v' - 4 A^^ cos 4 v') + — j 

BA,, BA,2 BA,3 ^ BA24I 

+ V, cos v' — — + 2 v, cos 2 v' — — + sin 3 v — — + sin 4 v -- — 

1 Ba ^ Ba Ba Ba J 


Bv BA Bv 

A + V (A , cos v' + 2 A cos 2v ) 

Ba Ba Ba •'* 
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^ (if - I §)] [-^1 + 2 ^ 0 ) + I q* sin ( 2 ^, + 2 ^„) 

"el + 4 v&o)j + i q 2 B-‘ |cos ( 2 v!-, + 2 ./.^,) ^ - | 2 ^j sin ( 2 ./-,+ 2 0 ,) 


f q^ cos (20^ + 20,) +^ q 2 cos (403 + 4^^) 


|q [sin (203 + 20,) - I sin (403 + 40,)j |j.| 



^^2 




A ( 20 s + 20 ^j) + ~ q 2 sin ( 2 v ^3 + 20 ^) 


64 


q 2 sin (4 0 ^ + 40 ^) 


- s' 
2 ^ 




+ (ai 







+ (- 2 aj )-‘/2 (^2 _<^|y/ 2 ^-iB-i 


avj 

SAj 

3vi 

aT ^''2 

3e 

^ 3e 


(A^j cos v' + 2 A 22 cos 2 v') 


- (A 2 j Vj sin V* 


+ 4 A 


22 Vj 


sin 2v' - 3 A 23 cos 3 v' - 4 A^^ cos 4 v') 



+ Vj cos v' + 2v cos 2v' — li + sin 3v* + sin 4v* 

oe 3e 3e 3 e 




0 J cos (203 + 20,) +1 q 2 sin (203 + 20,) 


■•p q^ sir* (405 + 40,) 


1 1 ^01 
4 ^ 2 ' (20s + 20 q) ~ 


201 sin (203 + 20 ,) 


- I q* cos (203 + 20,) + ^ q* cos (403 + 40,) ^ 




+ - q sin (203 + 20 ,) - sin (M + 40 ,) 


'o)l 

dej ’ 


I I 


I I 1 1 II 1 1 I II I I 


I I 


1 


B'Z '2 


= ~ J ®2‘ ( 2 >/.^ + 2 >/-o) + |; q 2 sin {2^p^ + 2 </-„) 


3 2 • 

sin 

64 




+ (- 2 “j) (»2 - ^-l g -1 


3 v 3 A 3 v 

Aj -r — +v - — + - — (A„ cos v' + 2 A„ cos 2 v') 
. ^-^0 


/ 3 Mg 3 v 

- (A«- V sin v' +4 A<,« v- sin 2 v* - 3 A„- cos 3 v' - 4 A., cos 4 v') - — + - — 

22 1 23 24 


3 A, 


+ V COS V* + 2 v, cos 2 v* ~ + sin 3 v* 

A /A'n 1 'A'n 




3^70 


3^0 




+ sin 4 v 


'!^1 

^^0 J 


3 t?o 


v//j cos ( 2 >/'^ + 2 i/>j,) 


+ -| sin ( 2 'Z'j + 20 (,) - ^ q 2 sin(4<p^ + 4 i/'p)j+ i q^ -jcos ( 2 '/'^ + 21^^) 


3 i//j 


| 2 '/>j sin ( 2 '/>^ + 2 >/>q) ' ( 2 '/’s + 2 '/’^) + ^ q^ cos ( 4 </<j + 


3 'Z's S'/'o 

'J VH " ^ 


3 

q 


sin (2'p^ -I- 2 si'o> - - sin (44i^ + 4 </>(,) 

O 


Bq 

3-^0 


Bi/f 

1 

3^0 

+ H 

3i/)j 

3^2 

Ba 

Ba 

Ba 

Ba 

Ba 

Bi/f 

B'/'s 


3^1 

3 v ^2 

Be 

Be 

Be 

^ Be 

Be 

Bi// 



3>Ai 

3 v ^2 



3’7 o 

3-^0 

3^0 


This completes the computation of the partial derivatives of the uniformizing variables E,v , 
and <P with respect to the orbital elements a , e , and % when the calculation is followed throv^h 
terms of the second order. If, however, second-order precision is not necessary, we can 
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eliminate the terms with the subscript ^^ 2 ^’ (thus omitting all partial derivatives of E29 ^2 y 
and 4^2)9 above partial derivatives of the uniformizing variables reduce to: 

Ba Ba «-Ba 


3e Be Be ’ 

BE ^ 38 3Ej 

^Vo 3% 3^0 

3L , !!li + !lo + ^ 

Ba Ba Ba Ba 

^ 

Be Be Be Be ’ 

Bv _ ^ ^ 

3 % ’ 

d<p ^ ^ ^ ^ 

Ba Ba Ba Ba 

^ ^ ^ 

Be Be Be Be ’ 

Bi/» _ 3i^ Si/'j 

H ^ H H ' 

We now continue with the necessary equations preparatory to the partial derivatives of the 
orthogonal coordinates X , Y , and z . 

^ ~ " ’To) x)** (1 - ’To sin^ , 

rlfl ” Thq 


^ = (1 - ’To) 'P (sin v)*‘ (1 - sin^ i/-)-s /2 ^ , 


(1 - '^o) sin i/<(sin y)'T (1 - sin2 ^)-3/2 


- i7o cos i/< sin* ■// (sin y)"‘ (1 - 75* sin* i/;)"*/* , 
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B <^ 


(1 - (1 - X + B 3 ^ vlvl"* sin 2 <l> (a| - a§)' 


3 a 

1/2 „ _ 3 .. 

3 a 




( ^ a 

'« r* 3 a ■ ^ 3 a/. 


+ ^ 3 .(°-i - ’Jo (1 - It 


Ba 


- ’ j ;^ (1 - (1 - ^ if ^ ^ ^ 3 , 


^’^2 „ 3 ii/ , ^^3 

-T-^ + B, + '/' ^~ 


^ 7 ^ 77 -s sin 2 s& ^ + -^ 7,2 r?-“ cos 24 > 


8 


3 a 


3 a, 
- 1/2 3 


Ba 1 

|(-2a,)-‘/^ ^ +a3(-2a^)-3/3 _J(A3V . A^ 

+ A 33 sin 3 v + Aj^ sin 4 v) - c^ S ^ “ 2 aj )"*^3 |^ 


sin V + A 32 sin 2 v 


(A 3 + cos V + 2 A 32 cos 2 v 


Bv ^^31 ^^3 2 

3 A,, cos 3 v -f 4 A-. cos 4 v) + sin v — — + sin 2 v 

33 oa 


Ba 


3 A „ 3 A3 

+ sin 3 v-g^ + sin 4 v — +v _ 


3 <^ 

3 e 


(1 - 7 ?^)''^^ (1 - X + B 3 /< +-^ 7)2 7 )-“ sin 2 si- 


( a2 _ a|)-‘/2 7 ,^ ^ 


- a 3 ( a 2 - a 2 )- V 2 


"• 3'''^2 3 


( ^“2 ^‘’ 3 \ 

’’O V“2 17 - -3 17jJ 


+ - “ 3 )' ’Jo 


( i -' il )-‘''< i - 7 -,=)-'’| f - 


- 7 ;- (!->«)-■'■ 0 - 7 ^)-" J , 


3 B, 


1 - 7)2 7,7 Sin 2 /^ ^ + A „2 7 )-“ cos 2 </< U 


Ba. Ba 

(- 2 a ,)-‘/2 _ ^ 2 a 3 >- 3/2 _ 


(A 3 V + Ajj sin V + A 32 sin 2 v 


+ A 33 sin 3 v + A 3 ^ sin 4 v) - c^ (- 2 aj )-*^2 


(A 3 + A 3 , cos V + 2 A 32 cos 2 v 
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+ 3 Ag 3 COS 3v + cos 4v) 


By 

le 


+ sin V 


3e 


3A,, 

+ Sin 2v — — + sin 3 V — + s in 4v 

de de 


Be 



d<t> 


(1 - 7t2)-i/ 2(1- 77-2)- 1/2 y +B 3 >/- + | 5 -^?-^i'*sin 2v/, 


“■3(°-2 - V, 


' 3 '^ ^32 
Ba^ Ba, 


(a2_c^V‘^^^o ^ 


■, 2 n- 1/2 


+ “3(^2 - 


L 0 '0 ^0 '0 


■■f ’’2® 2>A '^0 •^j'^cos 241 ^ + 773 (1 - ?72)-2/2 (1 _ ^-2)-i/2 


2 

+-YJ’ 7 o 2 ? 2 '' sin 2 ./- 


Ba_ Ba 

(-2a )-l/2 + a (_2a,)-2/2 1 

^ ^ 3-^oJ 


(A-V + A, j sin V 


-I- A 32 sin 2 V + A 32 sin 3v -h A 34 s in 4 v) - ( - 2a^) 


- 1/2 


(A 3 -t- A . cos V + 2 A 2 cos 2 v 


BV ^^31 ^^32 BA 33 BA 3 ^ 

+ 3A,,cos3v + 4 A- - cos 4v) + s in v — -h sin 2 v — + s i n 3v — -f s in 4v — — 

33 34 


^^3-1 


BX 

Ba 

BX 

Be 

BX 

^^0 


(p2 ^ c2)‘l ^ [ a e si n E + 1 _ e cos E] - (1 - Pq s i n^ i//) ^ 7)2 s i n i// cos i/; 

' Ba ‘ 


Ba 


- Y 


(p 2 + c2) 1 p a ( e s in E - cos E^ - (1 - s in^ v//) ^ 7^2 s in 0 cos \p 


B 0 

’ 

B0 

Be 


= X 


(p2 + c2) 1 p a e sin E 


be 


BY 

Ba 


= Y 


- (1 - sin 2 0)^1 t}q sin <// cos 0 

/ • IT 

( a e sin E — — 

\ Ba 


+ sin 4> 


- Y 


Bc^ 


(P + c2)-l I 


+ 1 - e cos E 


Be 


- (1 - p 2 sin 2 ,^^-1 7)2 sin 0 cos 0 

(p2 + c2)“^ pa sin E cos E^ 

- (1 - 7 ^ sin^ 4>)^^ Vq sin cos i/> 


+ X 


-h X 


Be/) 

Be 


51 




BE 


(fp‘ + c2)"^pae sin E 


~ (1 - 77^ sin2 i//)-i Tj^ sin \p ^77^ cos i// + sin j 


+ X 


B0 


BZ 

Ba 


?7 q (1 - e cos E) ^sin + a cos i/; ^ ^ ’ 


= 77^ (1 - e cos E) a cos \p ® ^ ’ 


a(l - e cos 
^^0 


E) ^sin ^ + Vq cos xp + a e 77^ sin i// sin E |^ . 


DIFFERENTIAL CORRECTION; TIME-VARYING PARTIAL DERIVATIVES 
WITH RESPECT TO ANGLE-EPOCH VARIABLES 

We now compute partial derivatives of the time -dependent parameters from the orbit 
generator with respect to the orbital elements /3j , and . This procedure is analogous to 
the one followed in the preceding section. Whenever a partial derivative with respect to j6^ is 
not given, it is assumed to be zero. 

= 2 m/ , 

d/s, 

^=~ “2* ®1 ®2‘ ’ 

dy«2 

^-/-s 


— 1 = 2771/ a"^ Aj* (a + bj + Aj) , 

d/3, 


dS . p.-i 

_=(!-. co.E) 


IL =(1 -e' cos8)-l ^ , 

B/S, d/ 3 ^ 


— 2. = (l-e^)sin£(sinv') ^(l-ecosS)~^ 

d/3, 


dfi, ' dfi. 
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= (1 - e2) sin 8 (sin v')'> (1 - e cos 8)"^ — , 

dyS^ 3/3j d/3^ 

^ A, b;' ^ . 


3 V 

— = (- 2 a ,)->/2 ( a ^ - a 2 ) l /2 ^-1 A, B'* ^ , 


d/3, 


dM, 

^=-(a.b,)-l 




(A, ^ 


2 


(-2a,)l/2 (a2 -a 2 )-l ''2 cos (2-/-, +2^,,) 



dM 


(Aj + c2t,2 A,B,B-‘) 


^^0 

3^, 


-•J ‘=^'^0 (- 20 - 1 )^^^ (a| - a2) 1/2 cos(2i/-^ + 20^) 




dE, 

^ = (1 - e' cos8)-i 
d/3. 


dMj 


[l-M^e'(l-e'cos8) ^sinS] 


- e' (1 - e' cos F) ^ 


Be 


[sin P 


+ (1 - e' cos P) ^ cos e - (1 - e' cos p) ^ sin^ p] , 


BEi 

_ = (1 - e' cos 8)*^ 

dp. 


dM, 


[1 - M^e'(l - e' cos p)"^ sin P] 




e' cos p)’^ 


BP 

d/3, 


sin P + ~ M, (1 
2 ^ 


e' cos P) ^ cos p 


— Ml e' (1-e' 


cos P)"^ sin^ P 
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(1 - e^) sin (S +Ej)(sin v")‘‘ 


[1 - e cos (S + Ej )]-2 


'Bg BEj 




= (1 - e^) sin(g + E.) (sin v")~* [1 

3^2 


e cos(g + Ej)]~* 


B 8 ^ 



B/ 3 , 


= (- 2 a,)- 


,-1/2 


^“2 


a2)t/2 



+ (A 21 COS v' + 2A^2 cos 2v') 


' SM3 avoV 


+ i B"‘ cos ( 2 <//^ + 2 '/'„) 


! V \ 

\ 3 / 3 , "b^,/’ 


B^, 

3^2 


^=(- 2 a ,)-‘/2 


("2 


a2)l/2 ^-lB-1 



+ (A2^ cos v' 


+ 2 A22 cos 2 v') 


/BM, ^ Bv^Y 
^ ^^ 2/- 


+ — q 2 Bj* cos ( 2 i/'j + 241^) 
4 


3 /S, B/ 3 , 


The following time- dependent partial derivatives with respect to the orbital elements , 
, and are computed only if the differential correction is carried through terms of second 
order: 


^^2 , r 

— 1= - (a + b.)“i <^A, — L . /A 


3 M- 3 v 

,, cos V' + 2 A„ cos 2 v') 


+ (_2a,)I/2 (a|_o.2)-l/2 


3'Z'l , 3.//, 


+ '/'i sin ( 2 iAj + 2i/<q) 


^ 

3 / 3 , 3 / 5,7 


_ iq 2 cos ( 2 i//, + 2 >/-g) 



3^1/ 


+ 


_1_ 

16 


q 2 cos + 40 q) 
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BM2 ^ f Bv /BM Bv.\ 

- (a . b,)- |a, — + (A„ cos V' . 2 A,, cos 2v') + — j 

+ c 2 t ?3 ( _ ^ 0,2 a 2 )-i /2 


Bi// 1 

^ (S'A, + 2v-„) ^ 


Bi/)j 


. sin (2^, . 2^„) y-. ^ ^ W. 


1 Z^’/' Bvi; 

+ — q2 cos (4^3 +40„) (^ 


3/?i 


= [1 - e' cos (g + Ej)] 


-1 —1 


3/^2 3/^2 

BM, 


3/3^ 


- 1^2 [1 - e' COS (8 + Ej)] ^ e' sin (8 + Ej) + 1 , 


Bg 


3/Si 3/3j 


3E, , BM, 

• — ? = [1 - e' cos (g + Ej)]-i — 


3A 


3A 


/Bg BE,\ 

-M [l-e'cos(8 +E-)]'^e' sin(8+E.) ( + •) , 

' 3 «y> 


[1-e 

' COS 

(S +Ei)] 

Bg 

BE^ 

3E, 

3/3/ 

3 / 3 , 

'3/3, ’ 

Bg 

BE, 

3E, 

3/3/ 

3 / 3 , 

' 3^2 ’ 


3^1 

3/3, 


- e2) sin E(sin v) ^ (1 - e cos E)“2 

BE 

3/3, ■ 

/BM^ 

y " 

BVq 3v,\ 

3/3, "3/i,; 


BE 

/BM, 

Bvq 3v,\ 

- e2) sin E(sin v)’^ (1-e cos E)"^ 

3^2 ■ 

\3^2" 

3/32 " 3 / 3 J 


Bv _ ^ 3 Vq 3v^ 3v^ 

^ B^ ^ ^ Bi3, 


Bv _ ^ ^ ^ 

^ ^ ^ ^2 '*’3/3/ 3^82 
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r Bv- 

(- 2a^)-‘/2 (a* - ^-ig-i |a, ^ ^ (A^j cos v' + 2A,j cos 2v') 


"(A, 


/ 3 M, 3 v„Y 

Vj sin v' + 4 Aj 2 Vj sin 2v' - 3 A^j cos 3v' - 4 A^^ cos 4v') 


+ -| ( 2'/'3 + 2%) - |2'/'j sin { 24 >^ + 2</-g) 


--i cos { 24 <^ + 20„) + — cos (4>//^ + 4'/'o) 


16 


^ 3 ;^' 

3/6, ^3^j 


3'/'2 


(-2a^rJ^2(a2-a2)i/2^-*B;‘ [a 2 (A 2 , cos v' + 2 A 22 cos 2v') 


3 / 5 , 


/ 3 M^ 3 vg VI 

- (A2, Vj sin v' +4 Aj2 Vj sin 2 v* _ 3 Ajj cos 3 v' - 4 A24 cos 4 v') + -g^jj 

+ lq2 B-i |cos (2^3 + 2^0) ^ - 1^2^, sin(203 + 2./-„) 


■-^ q2 cos (^2yp^ + 2</’o) q* cos + 4c/<j) 


3'Z's 3./<o 

+ 


3^2 ^^2 


3.// _ ^ ^ 

^ ^ ^ ^/ 5 i ^/ 5 i 

3,/, _ 3'/>3 ^ ^ ^ 

3^" 3/3j 3/6^ 

This completes the computation of the partial derivatives of the uniformizing variables 
E , V , and with respect to the orbital elements P2 > ^3 U^^ose with respect to p^ are all 

zero) when the calculation is followed through terms of second order. If, however, second- 
order precision is not necessary, we can eliminate the terms with the subscript ”2^' (thus 
omitting all partial derivatives of M2 , £3 , V2 , and ), and the above partial derivatives of the 
uniformizing variables reduce to: 

BE ^ Bg 

3/5j" 3yS, ^ 3/S, ’ 
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BE 

=11. 

BE 3 


BA 

BA 

BA ’ 


Bv 

3M, 

^''0 

Bvj 

BA 

"ba 

"ba 

"^A 

Bv 

BM^ 


Bv, 

BA 


" ^A 

ba 

3^ 

BA 

3A 

BA 

BA 

~BA 

^ BA 


B\/i 

BA 

BA 

BA 

BA 


BA 



We now continue with the necessary equations preparatory to the partial derivatives of 
the orthogonal coordinates X,Y, and z. 


(1 - rji) sin >/. (sin x)'* (1 - ^0 2^ , 

(1 - T 5 *) sin 4 > (sin x)"‘ (1 ~ ^ ’ 


B/ 3 j 


= -3 


(1 (1 -^2 ) 


2 x- 1/2 




+ B, 


^ + A. ^2 n-* cos 2-^ 

3 / 3 j 16 ® ^ 


'^xp 


- c* tt 3 (- (Aj + Ajj cos v + 2 Ajj cos 2 v 


+ 3 A 33 cos 3v + 4 A 3 ^ cos 4v) 


Bv 


BA ’ 




^0 


(1 - 


■1/2 


(1 - ^V) 


2\-l/2 


By 

BA 


BA 



^^p' 

7,2 7 ,-^ cos 2 ^ ^ 
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- c2ag(-2a^) (A3 + A3 j cos v + 2A32 cos 2v 


Bv 

3A33 cos 3 v + 4 A3^ cos 4v) 


, 

^ = X l(f^ + c 2 )-‘ p a e sin E 

d/?l L 

- (1 - vl sip 2 ./))-> sin 1/- cos i/< - Y , 


3 X 

3/3, 


X 


(p 2 + p a e sin E 


BE 

3^, 


- (1 - sin^ 


?7^ sin ip cos ip 




BX 

B/33 


- Y 


by 


= Y 


(p2 + c^) 1 


p a e sin E 


BE 

B/3^ 


- (1 - ■^0 sin .// cos i// ^ ^ ’ 

|^ = Y [(p^ + c2)-> paesinE|^ 

3 / 3 , L ‘^'“2 

- (1 - ^ sin2 172 sin ./-cos ./. + X , 


3 Y 

3/3, 


BZ 

3^ 


= a 77o j(l - 


e cos E) cos \p 


'dip 

W, 


+ e sin \fi sin E 
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BZ 

3/3, 


- a 77 


'o [d - 


e COS E) COS >// 


30 


+ e sin 0 sin E 


3E 

3A 



0 . 


THE EQUATIONS OF CONDITION 

Now that we have found 3x/Bq, , 3Y/3q. , and Bz/Bq. for the Izsak orbital elements q; (i = 1, 
2, . . . , 6), we can complete the differential correction process by determining the equations of 
condition. First we expand and substitute into the matrix relation given in ANALYTICAL PRO- 
CEDURE OF DIFFERENTIAL CORRECTION. The matrix relation, when expanded explicitly, 
yields the following eighteen equations: 






— — = - Sin 0 




+ cos 0 


sin« H 
° Ba 


+ COS 6^ 


BZ 

Ba 


Ba 


BX 

sin 0 cos 9 cos 0 

^ Ba ^ 


^ BY . . BZ 

COS^D ^ 


^^1 ; BX . , BY 

— = cos ^ + sin 0, — 


_=-sin>^x 


BX 

sin - — + cos 0^ sin 9 

"Be 


BY 

Be 


+ cos 9^ 


BZ 

Be 


Be 


= sin 0 cos 9^ 


BX . 

— - cos 0x 


° Be 


• a 3Z 




= cos 4i + sin 

3i7g BT7g 


3'^o 


-sin0 sin^j. 


BX 

3^0 


+ cos 0 sin (9 


BY 

3^0 


cos 


BZ 

3^0 


=:sin0 cos^^ 

3i7o 


- cos l// cos 

3-^0 


9 W « 

3-^0 


3\i 

— - = cos 

B/3j 


BX 


+ sin i//x 
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= - sin 0 sin 

3/3j 


— _ + COS 0^ Sin 
3/3j 





— = sin cos 0JJ 


3X 

3/3i 


- COS 0J^ COS (9jj 


BY 

3/3i 


+ sin 6 ^ 





cos 0^ 


BX 

^^2 


+ s in 0^ 


BY 


— ^ = -sin0sin^- 
B/32 ' 


BX 

^^2 


+ cos 0^ 


BY 

^^2 


+ cos 0 ^ 


BZ 


3Zm . , 

3 ^ = sin cos 6)jj 


3X 

3^2 




3Y 

3/^2 


+ s in 


BZ 


3/^3 


= cos 0^ 


BX 


+ s i n 0^ 


by 

3/^3 



sin0 sin 6^^ 

”x D 


BX 

^^3 


+ cos 0^ 


sin 


BY 

3^3 


^Zm . , . 3X 

5 ^ = 5 ^ 


- COS 0^ cos — 


by 

^^3 


The last two equations have only two terms on the right-hand side because of the fact that 
Bz/B /33 = 0. We can now write out explicitly the twelve coefficients to be inserted into the equa- 
tions of condition: 


3a 3X^, 3a 3Y^ Ba ^ 3Z„ 3a 

Be " 3X„ 3e " 3Y„ 3e " 3z„ 3e 

3L, _ 3L, 3X„ ^ 3L, 3Y„ ^ 3L, 3Z„ 

3t7q ■ 3X„ 3t7o 3Y„ Bt]^ ^ 3Z„ Bt,^ 

^ 

3/3j " 3X, B/3, " BY„ 3/3, " 3Z„ 3/3, ’ 

3/32 ■ 3X„ 3/3^ ■• 3Y„ 3/3, ^ 3Z„ B/3, ’ 
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BL, 3L, BX„ 3L, BY„ BL, BZ„ 

Ba BX„ Ba BY„ Ba BZ„ Ba ’ 

BM, BY„ BM, BZ„ 

Be BX„ Be '*’ BY„ Be BZ„ Be ’ 

BM, BM, BX„ BM, BY„ BM, BZ„ 

3’7o'^X„B77/by„B,?/bz„B7,3 ’ 

BM, BM, BX„ BM, BY„ BM, BZ„ 

3^1 -BX„B/3, "3Y„3^/3Z„3^i ’ 

BM, BM, 3X„ BM, 3Y„ BM, BZ„ 

B/?3 ' BX„B^3 "3Y„B/?3 "3Z,B/?3 ’ 

^ ^1£m 

9/^3 ■ B/33 "by„ 3^3 "bz„ 3/33 • 

Finally, the two equations of condition corresponding to each observation are given ex- 
plicitly by: 


BL BL^ 

AL = L„-L,=^Aa + — Ae+— — A/3,+ — A^ 3 + — A ^3 


3M 


3M„ 


BM, 


AM = M. - M, = — Aa + -A- Ae + ^ 

“ Ba Be drj^ 


BM^ 


BM, 

^2 


BM, 

3/^3 


FITTING BY METHOD OF LEAST SQUARES 

We have accumulated a set of 2n linear simultaneous equations in six "unknowns,” as 
follows: 


(AL). = (Lq), - (L = 


z: 

5 t= 1 



Aqj 


(AM), = (Mq), - (Mi), = 


E 



1, 2, • • • ,n 


where q.(j = 1, 2, . . . , 6) = a, e , /3^. We regard the Aq . as "unknowns," and the 

number, n, of observations in the set is fixed in advance (see above under THE STANDARD DEVIA- 
TION OF FIT). The above equations, written in matrix form, become 
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[Aqp = [AL.] 
[Aq^ = [AM^] 


where the matrices of partial derivatives have n rows and six columns, the matrices of un- 
knowns have six rows and one column, and the matrices of observational residuals have n 
rows and one column. Recall that, in the general case, the observed direction cosines (L^). 
and (^q). have associated weighting factors (wL)i and (w„)i, respectively (see above under 
COMPUTATION OF DIRECTION COSINES). 


For purposes of this section, it is unnecessary to distinguish between direction cosines 
L and M or between weighting factors wl and w„ . Further, it is not significant, for the present 
purpose, that the constant coefficients in the linear simultaneous equations have the form of 
partial derivatives. In order to simplify the notation in what follows, we combine the two 
matrix equations, each coefficient matrix having n rows, into a single matrix equation where 
the coefficient matrix has m = 2n rows. Then the matrix of constant terms (i.e., observational 
residuals) also has m rows. We rewrite the above two equations in the simple general form 


AX = B, 


where A = [a. .] has m rows and six columns and represents the coefficient matrix of partial 
derivatives, x = [x.] has six rows and one column and represents the matrix of unknowns, and 
B = [b.] has m rows and one column and represents the matrix of observational residuals. 

The number, m, of equations we obtain by expanding the matrix relation is generally much 
greater than the number (six) of unknowns, and since the observations contain inherent random 
and possibly systematic errors, no exact solution of the simultaneous set exists. According to 
the principle of least squares, the values of the unknowns x . which are preferred are those 
which cause the sum of the squares of the residuals after the fit to be a minimum. The so- 
called ^residuals after the fit" are calculated by substituting the approximate solution for the 
X. in the matrix X, and subtracting the matrix AX from B. When the equations of condition 
have different weights, the least- squares solution is that which minimizes the sum of the 
weighted squares of the residuals after the fit, where each square is multiplied by its corre- 
sponding weight. 

The least- squares criterion is satisfied by reducing the m equations of condition to six 
equations known as normal equations. This procedure, in which we adopt the usual notation for 
matrix elements (the first subscript denoting the row number and the second subscript the 
column number), is performed as follows. The first normal equation is obtained by multiplyii^ 
the first conditional equation by the second by the third by W 3 etc., and sum- 

ming the resulting m equations. The second normal equation is obtained by multiplying the 
first conditional equation by a^^, the second by third by etc., and summing 

the resulting m equations. If we repeat this process six times, we obtain the six normal equa- 
tions. It is seen that this process is equivalent to premultiplying the matrix equation AX = B 
by the weighted transpose of the matrix A, where the rows of the transpose are multiplied by 
the corresponding weighting factors. The set of normal equations can be represented by the 
new matrix relation CX = D , where 

m 

Cii = 2 ^ a^. (i, j = 1 . 2 . . 6 ) , 

k = l 
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and 


a,. b,w, (i = 1.2.---.6) . 

k = 1 

Of course, if the weighting factors are not present, i.e., = 1 for all k, the elements c. . 

are precisely those of A'^A and the elements d. are precisely those of a'*“b. Here the super- 
script indicates the transposed matrix. 

We now have a system of six equations in six unknowns, since C is a square (and sym- 
metric) matrix. In order to solve this system, we use a method known variously as the Gaus- 
sian elimination method or the method of pivotal condensation. This has the effect of reducing 
the square matrix to an upper triangular matrix (i.e., all elements below the principal diagonal 
are zero) which represents the same solution for the x . . 

To begin this process, we choose the element of the first column of matrix C greatest in 
absolute value, say . We then divide all the elements of the k^^ row (the ^^pivotal row^') by 
the so-called dominant element (or ”pivot”)? This done, we exchange the corresponding 

elements of the pivotal row with those of the first row. The leading element of the matrix, 
is now unity. We now replace all the elements in each row beginning with the second row by 
the following procedure: multiply all elements in the first (pivotal) row by the element in the 
first column of each row successively and subtract this product from the corresponding ele- 
ment of the successive rows. Mathematically, this is indicated by 

^ij = ^ij (i = 2, 3, • • *.6 ; j = 1 , 2, • • •, 6) . 

Since 1, it is obvious from this equation that = 0 for all i = 2, 3, . . . , 6. That is, all 
elements in the first column except for the first (diagonal) element are replaced by zeros. 
Essentially, we have added suitable multiples of the pivotal row to all the other rows so that in 
each resulting row the element in the first column vanishes. 

Consider the matrix with five rows and five columns obtained by deleting the first (pivotal) 
row and the first column. Now select as a new pivotal element the largest element in absolute 
value in the new first column of the five- by -five matrix and repeat the entire process with re- 
spect to the square matrix of order five. 

Continuing in this manner, we have finally a single nonzero (diagonal) element in the last 
row. The procedure is completed by dividing this final row by the diagonal element. The result 
is an upper triangular matrix with ones along the principal diagonal. Note that all operations 
described above to be performed on the original square matrix C are elementary row operations 
(i.e., an operation belonging to one of the three following types: the interchange of any two rows; 
the multiplication of a row by any nonzero constant; the addition of any multiple of one row to 
any other row). Thus, the triangular ization process does not change the solution to the simul- 
taneous set of linear equations as long as the operations performed on matrix C are performed 
in an analogous manner on the elements of the column matrix D . This can most readily be done 
by augmenting the six-by-six matrix C by a seventh column composed of the elements of D . In 
practice, the six-by-six matrix C is further augmented by a six-by-six identity matrix placed 
in columns eight through thirteen. The purpose of this is to determine ultimately the inverse of 
the coefficient matrix C, from which we may easily find the standard errors of the least- squares 
solution for the Aq .. Note that the various columns of C"^ can be found in succession by solving 
the matrix equation CX = l^ for the colunan matrix X , where I. represents the i column (1 = 1, 
2, . . . , 6) of the identity matrix of order six. We can thus view the six columns of the identity 
matrix placed in the augmented six-by- thirteen matrix as constant right-hand side column 
matrices replacing B in successive least-squares solutions of the matrix equation. These suc- 
cessive solutions are determined simultaneously in the Gaussian elimination method simply by 
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forming the augmented matrix E and performing the elementary row operations on all thirteen 
columns. The augmented matrix E appears as follows, after the normal equations are deter- 
mined, but before the elementary row operations are begun: 


n 

‘'la 

^13 

^14 

^IS 

= 16 

dl 

1 

0 

0 

0 

0 

0 

21 

22 

^ 23 

^^24 

'=25 

= 26 

^2 

0 

1 

0 

0 

0 

0 

31 

^32 

^33 

‘^34 

*= 35 

= 36 

^3 

0 

0 

1 

0 

0 

0 

41 

‘=42 

S3 

*^44 

^=45 

= 46 

^4 

0 

0 

0 

1 

0 

0 

51 

^52 

^ 53 

^^54 

*= 55 

= 56 

d5 

0 

0 

0 

0 

1 

0 

61 

^^32 

^ 63 

^64 

= 65 

= 66 

de 

0 

0 

0 

0 

0 

1 


After the triangularization process, the augmented matrix E is transformed to a matrix 
(call it F ) of the form: 


1 

fl2 

fl 3 

fl 4 

^15 

fl6 

^11 

gl2 

gl 3 

gl 4 

^15 

^16 

gl 7 

0 

1 

^23 

^24 

^25 

^26 

^21 

®22 

^23 

^24 

^25 

^26 

^2 7 

0 

0 

1 

^34 

^35 

^36 

^31 

^32 

^3 3 

^34 

^35 

^36 

^3 7 

0 

0 

0 

1 

^45 

^46 

641 

^4 2 

^43 

^44 

^45 

^46 

647 

0 

0 

0 

0 

1 

^56 

^51 

^52 

^53 

^54 

^55 

^56 

^57 

0 

0 

0 

0 

0 

1 

^61 

^62 

^63 

^64 

®65 

^66 

^67 


The first six columns of F represent the triangularized coefficient matrix and the re- 
maining seven columns represent successive constant right-hand side matrices, each of which 
is associated with a particular column solution matrix. At this point, it is only natural that we 
augment the column solution matrix X (corresponding to the seventh column of F only) to a 
six- by- seven solution matrix which contains X as its first column. The remaining columns 
of Y will contain the inverse of the coefficient matrix C of the normal equations. 

We can now write explicitly the set of linear simultaneous equations in the triangularized 

form. 


Vu + ^12 ^21 + ^3 y3i + ^4 ^41 + ^5 
y 2 i + ^23 y 3 i + ^24 y 4 i + ^25 
y 3 i + ^34 y 4 i + ^35 
y4i + ^45 


ysi + ^6 yei 
ysi + ^26 yei 
ysi + ^36 yei 

ysi + ^46 yei 
ysi + fee yei 
yei 


gji . 

^ 3 i ’ 

^ 4 i ’ 

^ 5 i ’ 

^ 6 i 
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In the above, the subscript ^’i" assumes values from one to seven, corresponding to various so- 
lutions for the seven right-hand side sets of constants. 

The construction of this triangular system of equations is known as the forward solution 
and the process of obtaining its solution is called back- substitution. The last equation in the 
triangular system gives the value for y^. directly. If we insert its value in the previous equa- 
tion, we can obtain , and so on for the remainder of the unknowns. Mathematically, the re- 
lation is 


and 

6 

Vji = gji - 22 ’ 

k = j + 1 

where j == 5, 4, 3, 2, 1 (in that order) and i = 1, 2, . . . , 7 (in any order). 

We have now completed the determination of the Aq. = y.^ by the least-squares principle. 
Theoretically, this procedure may always be followed to a successful conclusion provided that 
the m original equations of condition are independent; that is, provided that the determinant of 
the coefficient matrix c does not vanish. 

The formal solution by the method of least squares is now concluded, but ordinarily a 
measure of the ''goodness" of the least- squares fit is desirable. The residuals after the fit 
are assembled in the so-called residual matrix U, equal to B - AX. In terms of elements, 

6 

u. = b. - a. . X. (i = 1, 2, • ■ *,m) . 

j = 1 

From this, it is obvious that the sum of the squares (unweighted) of the residuals after the fit 
is given by 



We can now easily find the so-called variance -covariance matrix of the fit from the in- 
verse of the coefficient matrix in the normal equations. Recall thatC"^ occupies columns 
two through seven of matrix Y. The variance-covariance matrix is obtained simply by multi- 
plying each element in by the sum of the squares of the residuals after the fit and dividing 
this product by m-6 (the excess of simultaneous equations of condition over the number of in- 
dependent unknowns). If we represent the variance-covariance matrix by V, then we have 


where 






(i, j = 1,2, ,6) 



; = t 



By comparison with computations performed above in THE STANDARD DEVIATION OF FIT, 
we can see that the quantity is a standard deviation of fit. More precisely, is the 
standard deviation of the residuals after the fit, or the standard deviation of the least squares 
fit. It is not to be confused with 
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in the earlier notation, which is the standard deviation of the observational residuals, or the 
standard deviation of the observational fit. 


Finally, we can find the so-called standard errors, m-, of the six unknowns, Aq. =y.j. 
These are simply equal to the square roots of the diagonal elements in the variance-covariance 
matrix, or 

(i = 1 - 2 .- •• -, 6 ) , 

where Yj , j +i is the term on the principal diagonal of the inverse of the matrix C , corre- 
sponding to the unknown x . = y. ^ . 

ITERATIVE LEAST-SQUARES PROCEDURE 

The procedure for producing a mean set of orbital elements is essentially an iterated 
least-squares fitting of the first-order Taylor differential expansion of the conditional equa- 
tions to numerous observational values. Using the values for the Aq . determined by the method 
of least squares, as described in the preceding section, we can calculate the corrected Izsak 
orbital elements. 


a' = a + Aq^ = a -i- Aa , 
e' = e + Aq^ = e + Ae , 

^0 = ”^0 = ^0 * 

== /^i + Aq^ = /^i + ’ 

^2 =^2 =^^2 +^^2 ’ 

^3 ~ ^ 3 ~ ^3 3 ' 

At this point, it is useful to check that the improved or corrected elements are physically mean- 
ingful. For instance, it should be ascertained that the semi- major axis a' > 1 earth radius, that 
the eccentricity e' > 0, that the sine of the inclination, 77^, is not greater than unity in absolute 
value, and so on. 

It is now necessary to update the other parameters used in the differential correction 
process, based upon the improved orbital elements. Accordingly, the various parameters in- 
cluded under PRIME CONSTANTS II are re-evaluated using the improved set of elements. This 
done, the various parameters included under MUTUAL CONSTANTS are similarly re-evaluated. 
Now, assuming that the times of the various observations in the data deck are availably as 
needed, the Orbit Generator may be used to produce the required calculated values of the 
position and velocity components. From these components, we calculate the local coordinates 
of the satellite and then the computed values of the direction cosines (refer to COMPUTA- 
TION OF DIRECTION COSINES). Finally, the observational residuals are calculated. Thus, 
there are associated with each observation time in the data deck the following correspondir^ 
principal time-dependent quantities: 

(a) position and velocity components X,Y,Z,x,Y,Z; 

(b) local coordinates of the satellite x^,Yj^ ,z^ ; 

(c) computed values of the direction cosines ,M^ ; and 

(d) observational residuals AL,AM . 
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Next, the statistical analysis is repeated (refer to THE STANDARD DEVIATION OF 
FIT) wherein the following quantities are determined: the average observational residual, 
the standard deviation of the observational residuals from their mean value, the standard 
deviation of (the observational) fit, the upper and lower range limits for the observational 
residuals, and the standard deviation of fit of the accepted observational residuals. Once these 
quantities are found, the differential correction may be repeated. Of course, the time-inde- 
pendent partial derivatives are computed once only, while the time-varying partial derivatives, 
both with respect to the energy- momenta variables and to the angle-epoch variables, are com- 
puted for each observation time in the data deck. A new set of equations of condition can then 
be assembled and the fitting by least squares repeated. 

In summary, the following sequence of steps represents the iterative least- squares pro- 
cedure in producing a mean set of orbital elements for a given time span represented by a set 
of observation points: 

1. Correct the six Izsak orbital elements utilizing the values determined by the method 
of fitting the equations of condition by least squares. 

2. Update the parameters included in PRIME CONSTANTS EE. 

3. Update the parameters included in MUTUAL CONSTANTS. 

4. Produce sets of position and velocity components for each observation time using the 
Orbit Generator. 

5. Calculate the local coordinates of the satellite at each observation time. 

6. Compute the direction cosines of the satellite at each observation time. 

7. Determine the observational residuals for each time point. 

8. Perform a statistical analysis of the observational residuals to find various standard 
deviations and a statistically valid range within which observational residuals mnct fall for in- 
clusion in the fitting process. 

9. Begin the differential correction process by evaluating the time-independent partial 
derivatives. Then evaluate the time- varying partial derivatives for each observation point, 

10. Assemble the set of equations of condition. 

11. Fit the equations of condition by the method of least squares. First determine the six 
normal equations, then triangularize the system by the Gaussian elimination method, and finally 
use the back- substitution method to find the solution. 

12. Measure the ’’goodness” of the least-squares fitting by finding the residuals after the 
fit, the variance- covariance matrix, and the standard errors of the unknowns. Return to step 
number one. 

DEFINITIVE ORBITAL PARAMETERS 

The iterative least- squares procedure is generally terminated in one of two ways. Either 
the total number of iterations through the least squares routine is prescribed in advance, or the 
standard deviation of the observational fit is used as the criterion in halting the iterative method. 
If this standard deviation falls below a value prescribed in advance during a given iteration, 
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the precision of the differential correction is deemed sufficient at that point. Of course, both 
methods of terminating computation can be used concurrently; i.e., if the standard deviation 
does not meet the prescribed criterion by the p-th iteration, the differential correction 
process is halted. 

At the conclusion of the differential correction, the following definitive orbital param- 
eters may be found: 

The semi- major axis is found by multiplying a by the proper length conversion constant 
(3963.339 miles per earth radius or 6378.388 kilometers per earth radius). 

The eccentricity of the orbit is given by e. 

The inclination of the orbital plane to the equator is given by arcsin (0° ^ < 180°). 

The time of passage through the perigee point is found by multiplying by the proper time 
conversion constant (13.4472 minutes per Vanguard unit of time or 806.832 seconds per Vanguard 
unit of time). The time of perigee passage is given with respect to the reference (or epoch) time 
t 0 , which is that used as a basis for the observational times and that corresponding to the initial 
position and velocity components, X. ,y. ,Z.,X.,Y., Z.. 

The argument of perigee (measured in the orbit plane from the node to the perigee point) 
is found by multiplying by the angular conversion constant 57.295780 degrees per radian. 

The right ascension (measured in the equatorial plane from the vernal equinox) of the 
ascending node is found by multiplying /3^ by the angular conversion constant. (Note that these 
last two parameters are angles usually given as greater than or equal to 0° and less than 360°, 
so that some multiple of 360° may have to be added or subtracted to bring the values into this 
principal range.) 

The height of the perigee point above the earth^s surface is found by multiplying a(l-e) - 1 
by one of the length conversion constants given above. 

The height of the apogee point above the earth^s surface is found by multiplying a(l+e)-l 
by erne of the length conversion constants. 

The anomalistic mean motion is found by multiplying ^ angular conversion 

constant and dividing by one of the time conversion constants (this gives the mean motion in 
degrees per minute or degrees per second). 

The anomalistic period is found by multiplying 2 7ra^ 2 ^y time conversion 

constants. 

The mean anomaly (at the time of perigee passage) is found by multiplying - a " ^ 2 ^y 
the angular conversion constant. This expression assumes that the reference (epoch) time 
is zero; in general, the mean anomaly is found by multiplying -a ~^ '2 + t^) by the angular 

conversion constant. 


RESULTS OF PRELIMINARY APPLICATIONS 

Both the orbit generator portion and the differential correction process by least- squares 
fitting have been tested independently by application to actual satellite orbits. Primarily, use 
has been made of two relatively close-in yet drag-free satellite orbits, so that neither atmos- 
pheric drag nor luni- solar perturbing forces would exert major disturbing influences. The 
ANNA IB satellite (international designation 1962 BM 1; NASA identification number 56017) 
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was launched in October 1962 under the project direction of the U. S. Navy from the Atlantic 
Missile Range into a near-circular orbit of medium inclination. Its purpose was predominantly 
that of geodetic investigation. The Relay 2 satellite (international designation 1964 3A; NASA 
identification number 64031) was launched in January 1964 under the project direction of the 
National Aeronautics and Space Administration from the renamed Eastern Test Rar^e into a 
relatively high-eccentricity orbit. Its function was that of active -repeater communications 
satellite. Initial orbital parameters for Doth these satellites are given in Table 1. The obser- 
vational data for the Relay 2 satellite consist of direction- cosine pairs reported from fifteen 
tracking stations in the Minitrack network operated by NASA, while the data for ANNA 
IB consist of right ascensions and declinations reported from twelve stations in the optical 
camera network operated by the Smithsonian Astrophysical Observatory. It might be noted that 
no weighting factors were associated with any of the sets of observational data for either the 
Relay 2 or the ANNA IB satellite in the applications described in this section. 


Table 1 

Initial Orbital Parameters for Satellites Used 
in Preliminary Applications. 


Orbital Parameter 

ANNA IB 

Relay 2 

Perigee (statute miles) 

670 

1298 

Apogee (statute miles) 

728 

4606 

Period (minutes) 

107.8 

194.7 

Inclination to earth^s 
Equator (degrees) 

50.1 

46.0 

Semi-major axis (units of 
Earth's equatorial radius) 

1.1764 

1.7448 

Eccentricity 

0.00622 

0.23918 

Sine of the inclination 

i 

0.7672 

0.7193 


In order to gauge the intrinsic accuracy of the orbit generator, a double- precision ninth- 
order Cowell numerical integration program was utilized. Two numerically integrated com- 
parison ephemerides were produced, one using recently determined geodetic values for 
the zonal harmonic coefficients in the expansion of the geopotential and the other, the 
corresponding values for these coefficients based upon the Vinti potential (Table 2). The 


Table 2 

Zonal Harmonic Coefficients in the Geopotential 
Function Used in Generation of Numerically 
Integrated Comparison Ephemerides. 


Coefficient 

Geodetic Value 

Vinti Potential Value* 

j, 

1 . 0823x10’^ 

1.0823x10-3 

J3 

-2.3x10-® 

0 

J4 

-1.8X10-6 

-1.2X10-® 

J (n > 5) 

<1x10-6 

<lxl0-« 


♦The zonal harmonic coefficients for the Vinti potential function arc obtained from the 
relations: J ~ (-1)*^+^ and J 2 U +1 6* 
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numerically integrated ephemeris produced by 
the geodetic values of the zonal harmonic coef- 
ficients was used as a basis for comparison 
with both the numerically integrated ephemeris 
produced by Vinti values of the zonal harmonic 
coefficients and the ephemeris produced by the 
orbit generator based upon the Vinti potential 
function. Figure 1(a) illustrates the residuals 
of the X-coordinate between (1) the Vinti ephem- 
eris and the numerically integrated ephemeris 
produced by geodetic values, and (2) the numeri- 
cally integrated ephemeris using the Vinti values 
and the numerically integrated ephemeris produced 
by geodetic values. Figures 1(b) and 1(c) do like- 
wise for the residuals of the Y-coordinate and 
Z-coordinate, respectively. 


VINTI EPHEMERIS 
NUMERICALLY INTEGRATED 
EPHEMERIS USING VINTI 
COEFFICIENTS 


3 ^ \ i \ comparisons illustrated in Figure 1 

^ "I 0 \\ y based on the implicit assumption that the 

_i ^ initial position and velocity conditions do not 

3 I vlx contain any inaccuracies. In actual practice, 

^ - 20 - such inaccuracies are always present, and they 

(b) must be removed by utilizing observational data 
^ ^ ^ ^ in the differential correction. Figure 2 illustrates 

the determination of a mean set of Izsak orbital 

r fphfmfri<; elements by an iterated least-squares fitting of 

KiffAiiPoirAnv iMTFrpATPn differential solution to observational data for 

^ ^ 20 “ EPHEMERIS USING VINTI the ANNA IB satellite. In all cases, the total 

5 ? COEFFICIENTS / number of iterations through the least squares 

o ^ /> \ / fitting routine is prescribed in advance to be 

8 ^ Q \ This number is sufficient to attain conver- 

^ \ / \ / gence within a very small tolerance. In the 

^ - 10 - \y \ / ( X graph of each of the six orbital elements, three 

2 \ "curves'^ (actually a sequence of connected line 

a ^ segments) are shown, corresponding to various 

“ _ 3 Q _ numbers of observations included in the fitting. 

! ! ! ! An equation of condition results, of course, from 

0 100 200 300 400 500 6 mi -observation": either a single right 

TIME (minutes) ascension or a single declination value. One 

curve represents the minimum number of 

^ ^ , . . 1 r V#. . .1 equations of condition for a true least- squares 

Figure 1 -Coordinate residue s for Vinti potential fitting, viz., seven. This is associated with an 
ephemeris and for numerically integrated ephemeris observational arc length of approximately 45 
using zonal harmonic coefficients of the Vinti potential ^ represents twenty equa- 

each compared with numerically integrated ephemeris condition, or an addition of eleven equa- 

produced by geodetic values. tions, extending the observational arc length to 

approximately 75 hours. The third curve repre- 
sents fifty equations of condition, or a further 
addition of thirty equations, extending the observational arc length to a total of approximately 98 
hours. The starting point of each of the three arcs is the same, so that they overlap in time. Notice 
that each observational arc produces a somewhat different set of mean orbital elements, depending 
upon the additional observational values introduced. Physically, this may be explained as the re- 
sultant effect of forces not accomited for in the anal 5 d:ical theory. For example, electromagnetic 
disturbances, solar radiation pressures, aerodynamic drag, meteoric bombardment, etc., all influ- 
ence the mean set of orbital elements to the extent that they are reflected in the observational values. 
In performing the iterated least-square fittings, all the residuals corresponding to the preselected 


TIME (minutes) 

Figure 1 —Coordinate residuals for Vinti potential 
ephemeris and for numerically integrated ephemeris 
using zonal harmonic coefficients of the Vinti potential 
each compared with numerically integrated ephemeris 
produced by geodetic values. 
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SEVEN EQUATIONS OF CONDITION (45 HR ARC) 

TWENTY EQUATIONS OF CONDITION (75 HR ARC) 

fifty EQUATIONS OF CONDITION (98 HR ARC) 

Figure 2— Convergence of a mean set of Izsak orbital elements in ten iterated least-squares fittings of the differen- 
tial solution to observational data for the ANNA 1 B satellite, for various numbers of equations of condition. 


observation times were accepted at each fitting. That is, the acceptable range of values for the 
residuals constituted infinitely wide bands on either side of the mean value of the residuals. 
Mathematically, using symbols introduced in THE STANDARD DEVIATION OF FIT, 


[rj, r^] = 1 im [R - j cr, R + j cr] . 
j-oo 

Figure 3 illustrates the determination of a mean set of Izsak orbital elements by an iter- 
ated least- squares fitting of the differential solution to observational data for the Relay 2 satel- 
lite. In this case, however, the observational arc length and the total number of observations 
are held fixed, while the acceptability criterion for the observational residuals is varied. The 
arc length in all cases is one week, representing a total of eighty observations or a maximum 
of 160 possible conditional equations. In each of the six graphs, one curve corresponds to a 
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TIME OF PERIGEE PASSAGE, /3 (Vanguard units of time) SEMI -MAJOR AXIS, a, (earth radii) 








ONE-SIGMA RESIDUAL ACCEPTABILITY CRITERION 

TWO-SIGMA RESIDUAL ACCEPTABILITY CRITERION 

THREE- SIGMA RESIDUAL ACCEPTABILITY CRITERION 


Figure 3— Convergence of a mean set of Izsak orbital elements in ten iterated I east -squares fittings of the differen- 
tial solution to observational data for the Relay 2 satellite^ for various observational residual acceptability half- 
width criteria. 
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’^three- sigma” criterion, i.e., 


r^] = [R - 3cr, R + 3ct] , 

where, of course, a is the standard deviation of the observational residuals from their mean 
value. A second curve corresponds to a two-sigma criterion, and the third curve to a one- 
sigma criterion. Each curve is terminated when convergence of the orbital element is attained. 
Notice that convergence appears to be a slower process with a one- sigma criterion than with 
either a two- sigma or a three- sigma criterion. The rate of convergence in these latter two 
cases seems generally about the same. One effect of a wider acceptance range appears to be 
greater fluctuations in the value of an orbital element early in the iterated fitting procedure, 
although this is not always true. Also, despite the fact that differing numbers of conditional 
equations are accepted in the fittings depending upon the criterion for the residuals, the values 
of the orbital elements at convergence are remarkably similar. Refer to Table 3 for precise 
values, including the required number of iterations to attain convergence in each case. The 
uncertainties in the final significant figures (stated as ±X) are estimates based upon slight 
fluctuations in the values of the orbital elements in least- square fittings after convergence is 
attained. 


Table 3 

Values at Convergence of Izsak Orbital Elements for Varying 
Observational Residual Acceptability Half-Widths. 


Orbital Elements* 

One-Sigma 

Criterion 

Two -Sigma 
Criterion 

Three-Sigma 

Criterion 

Semi-major axis, a 

1.7444277 
± 1 

(12) 

1.7444278 
± 0 

(7) 

1.7444278 
± 0 

(7) 

Eccentricity, e 

0.2391624 
± 4 

(21) 

0.2391728 
± 3 

(7) 

0.2391818 
± 2 

(7) 

Sine of inclination, 

0.7231110 

i 2 

(21) 

0.7231020 
± 4 

(7) 

0.7231015 
± 2 

(9) 

Time of perigee passage, 

0.099885 

i 6 

(22) 

0.099985 
f 7 

(9) 

0.099942 
± 9 

(8) 

Argument of perigee, 3^ 

3.222265 
i 5 

(21) 

3.222232 

± 5 

(9) 

3.222278 
± 6 

(7) 

Right ascension of 
ascending node, 

‘ -2.380274 
± 3 

1 

(22) 

-2.380244 

• ^ 

(8) 

-2.380241 
t 4 

1 

(8) 


*Units of all elements are canonical (a in earth equatorial radii; /3 ^ in Vanguard units of time; ,3^ and in radians). 
The integers in parentheses refer to the number of iterations required to attain the converged value given. 


Table 4 presents the same information relative to the orbital elements as Table 3, but 
for an acceptance criterion fixed at two sigma, with the maximum possible number of condi- 
tional equations varied. The observational arc length remains one week, but the 160 maximum 
possible number of conditional equations are first reduced to one hundred, and then this num- 
ber is in turn reduced to forty. An attempt was made to maintain an even distribution of the 
observations throughout the seven-day period, while still operating on the "subset principle” 
(i.e., the set of twenty observations is a subset of the set of fifty observations, which is in turn 
a subset of the original set of eighty observations). 

Table 5 again records the same information relating to the orbital elements, but this time 
the parameter involving the order of precision in the differential correction is varied. Here the 
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Table 4 

Values at Convergence of Izsak Orbital Elements for Varying 
Numbers of Observational Points. 


Orbital Elements* 

40 Conditional 
Equations 

100 Conditional 
Equations 

160 Conditional 
Equations 

Semi-major axis, a 

1.7444276 

1.7444279 
± O'' 

1.7444278 
± 0 ^ ^ 

Eccentricity, e 

0.2391425 
± 5 ^ 

0.2391736 
± 2 

0.2391728 
± 3 ' ^ 

Sine of inclination, 77 ^ 

0.7231009 
± 3 

0.7230975 
± 3 

. 0.7231020 
± 4 

Time of perigee passage, 

0.100116 
± 9 

0.100002 , . 
± 7 

0.099985 
± 7 

Argument of perigee, 

3.222097 
± 5 ' ^ 

3.222261 .Qv 
± 6 

3.222232 
± 5 

Right ascension of 
ascending node, p^ 

-2.380249 
± 3 

-2.380248 , . 

± 3 

-2.380244 
± 4 


*Units of all elements are canonical (a in earth equatorial radii; /3^ in Vanguard units of time; and in radians). 
The integers in parentheses refer to the number of iterations required to attain the converged value given. 


Table 5 

Values at Convergence of Izsak Orbital Elements for Varying 
Orders of Precision in Differential Correction. 


Orbital Elements* 

First Order 

Second Order 

Semi- major axis, a 

1.7444276 , . 

± 0 

1.7444276 

Eccentricity, e 

0.2391427 - . 

± 8 ^ ’ 

0.2391425 

Sine of inclination, 

0.7231010 
± 3 

0.7231009 
± 3 

Time of perigee passage, 

0.100115 

± 8 

0.100116 
± 9 ^ 

Argument of perigee, P 2 

3.222098 , , 

± 6 ^ 

3.222097 . . 

± 5 ^ ^ 

Right ascension of 
ascending node, p^ 

-2.380250 ,- 2 ) 

± 4 ^ ' 

-2.380249 
± 3 


*Units of all elements are canonical (a in earth equatorial radii; /3^ in Vanguard units of time; j32 
and /3^ in radians). The integers in parentheses refer to the number of iterations required to attain 
the converged value given* 
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maximum possible number of conditional equations covering the one-week observational arc is 
held constant at forty, and the acceptance criterion is fixed at two sigma. The inexact designa- 
tions 'Tirst order^’ and "second order" indicate whether or not terms of purely second order 
are retained in the differential correction (refer to ANALYTICAL PROCEDURE OF DIFFEREN- 
TIAL CORRECTION). It is seen that retaining terms of purely second order adds immeasurably 
to the precision of the final converged results in all cases, and, similarly, does not affect the rate 
of least- squares convergence. 

The remainder of the figures display the convergence of perhaps the most significant 
single parameter in evaluating the efficacy of the differential correction process, viz., the 
standard deviation of fit. Actually, there are two standard deviations shown in each graph. 

The upper curve corresponds to a standard deviation of fit which includes all of the observa- 
tional residuals, while the lower curve corresponds to a standard deviation of fit which in- 
cludes only the observational residuals accepted at each fitting. Plotted on the same abscissa 
is a curve showing the number of equations of condition (or, equivalently, the number of ob- 
servational residuals) accepted at each iteration of the fittir^ process. 

Figures 4, 5, and 6 illustrate, respectively: the standard deviations for a maximum of 40, 
100, and 160 possible conditional equations covering the same one-week observational arc for the 
Relay 2 satellite. Note that convergence using a two-sigma criterion for the residuals, as shown 
in Figures 4(b), 5(b), and 6(b) is much more rapid than the convergence using a one-sigma cri- 
terion shown in Figures 4(a), 5(a), and 6(a). However, the convergence is not always so smoothly 
monotonic in the case of the wider acceptance range. Both these facts confirm what was said 
earlier about the convergence of the orbital elements. 

Table 6 supplies the values of the standard deviations at convergence for the various runs 
illustrated in Figures 4, 5, and 6, as well as several others not graphed. It also gives the number 
of accepted residuals at convergence and the number of iterated fittings required to achieve 
convergence in each case. 
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Figure 4 —Standard deviations of the observational 
residuals and the number of equations of condition ac- 
cepted at each iteration of the fitting process for a 
maximum of 40 possible conditional equations covering 
a one-week observational arc for the Relay 2 satellite. 


Figure 5— Standard deviations of the observational re- 
siduals and the number of equations of condition ac- 
cepted at each iteration of the fitting process for a 
maximum of 100 possible conditional equations covering 
a one-week observational arc for the Relay 2 satellite. 
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Figure 6—Standard devioHons of the observational re- 
siduals and the number of equations of condition ac- 
cepted at each iteration of the fitting process for a 
maximum of 160 possible conditional equations covering 
a one-week observational arc for the Relay 2 satellite. 


Figure 7 illustrates the standard deviations 
for a maximum of one himdred possible conditional 
equations covering an observational arc of only 
three hours for Relay 2. This is the three -hour 
period immediately following insertion of the 
satellite into orbit, when observations are re- 
corded at very frequent intervals in order to 
insure a wealth of data for the real-time differen- 
tial correction. Here, using a one-sigma criterion, 
convergence of the orbital elements occurs after 
only four (in some cases, five) iterations. The 
standard deviations of fit converge after three 
iterations to values of 0.425x1 0'^ (all one hundred 
observational residuals) and 0.145x10 (includ- 
ing seventy-seven accepted observational resid- 
uals). The graph shows that the standard devia- 
tions remain essentially constant after the third 
iteration and this is confirmed by the insignificant 
fluctuations in the orbital elements after the third 
iteration, although a total of ten iterations through 
the least squares fitting routine was prescribed 
in advance. 

Figure 8 illustrates the standard devia- 
tions for a maximum of one hundred possible 


Table 6 

Values of Standard Deviations of Fit and Number of Accepted 
Conditional Equations at Convergence for Various Runs 
Covering a One-Week Observational Arc. 


Description of Run 

Standard 
Deviation 
of Fit 
(all) 

Standard 
Deviation 
of Fit 
(accepted) 

Accepted 

Residuals 

Percentage 
of Total 

Iter- 

ations 

Total 

Conditional 

Equations 

Residual 

Criterion 

Order 
of D.C. 

40 

la 

2nd 

0.432 

0.157 

37 

92.5 

26 

40 

2a 

2nd 

0.440 

0.177 

38 

95 

12 

40 

2a 

1st 

0.438 

0.175 

38 

95 

12 

100 

1 a 

2nd 

0.373 

0.160 

95 

95 

11 

100 

2a 

2nd 

0.377 

0.169 

97 

97 

9 

160 

1 a 

1st 

0.307 

0.132 

141 

88.1 

24 

160 

la 

2nd 

0.307 

0.132 

140 

87.5 

22 

160 

2a 

2nd 

0.313 

0.172 

157 

98.1 

8 

160 

3a 

2nd 

0.315 

0.187 

158 

98.75 

9 


Note: All standard deviations of fit arc given in mils (i.e., in units of 10“^). The parenthetical word '*all” signifies that 

all of the observational residuals, AL and A M, were included in determining the standard deviation of fit; "accepted* 
means that only the observational residuals corresponding to the accepted conditional equations were included in 
determining the standard deviation of fit. 
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Figure 8— Standard deviations of the observational re- 
siduals and the number of equations of condition ac- 
cepted at each iteration of the fitting process for a 
maximum of 100 possible conditional equations cover- 
ing two distinct observational arcs for the ANNA 1 B 
satell iie. 


Figure 7— Standard deviations of the ob- 
servational residuals and the number of 
equations of condition accepted at each 
iteration of the fitting process for a max- 
imum of 100 possible conditional equations 
covering a three-hour observational arc 
for the Relay 2 satellite. 


cepted observational residuals) for Figure 
the values are 6.684 milliradians (all one 
radians (including ninety-five accepted observational residuals). 


conditional equations for two distinct nonoverlapping 
observational arcs for the ANNA IB satellite. Figure 
8(a) covers an arc of approximately nine days and 
fifteen hours, while Figure 8(b) covers an arc of 
approximately six days and nineteen hours. Both use 
a one-sigma criterion, and convergence of the stand- 
ard deviations occurs after five iterations in both 
cases. The values are a relatively large 38.379 
milliradians (all one hundred observational residuals) 
and 10.919 milliradians (including eighty-eight ac- 
8(a). For the somewhat shorter arc in Figure 8(b), 
hundred observational residuals) and 0.726 milli- 


The totality of data presented herein represents a small sampling of the preliminary 
applications by which the orbit generator and differential correction have been tested. Yet 
this sampling is indicative of the utility of the spheroidal method for artificial satellite orbits. 


CONCLUDING REMARKS 

The method of solution for unretarded satellite orbits discussed in this paper has been 
programmed, primarily in the FORTRAN language, for use on the IBM 7094 digital electronic 
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computer. It requires a relatively small number of computer storage locations, and the ana- 
lytical nature of the entire procedure assures a very rapid computational process. Extensive 
tests have indicated a capacity for generating coordinate and velocity points, based upon a set 
of empirically estimated initial conditions, in impressively short intervals of computer oper- 
ating time. 

Presently, work is underway on slightly modifying the accurate reference orbit to account 
for the effects of the most important perturbations of the neglected zonal harmonics, notably 
the third and the residual fourth. The inclusion of these perturbative effects by a procedure 
described in a recent paper by Vinti* is expected to improve the accuracy of the method so as 
to provide computed values agreeing with observation over a longer interval of time. 

In the future, a method of modifying the spheroidal potential for an oblate planet in order 
to permit the exact inclusion of the effects of the third zonal harmonic in the reference orbit is 
anticipated. Preliminary investigations are also being conducted into accounting for the luni- 
solar forces and aerodynamic drag. Further results will be published as they become available. 
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Appendix A 

Modifications for Use of Right Ascension-Declination Data 


Herein we present the modifications which must be introduced in order to utilize an alternate 
form of satellite tracking data known as right ascension-declination data. Such data are recorded, 
for instance, by the optical Baker- Nunn cameras of the Astrophysical Observatory of the Smithsonian 
Institution. In utilizing this data, the modifications to be described replace the material presented 
in the main body of this report imder COMPUTATION OF DIRECTION COSINES. 

A set of observation data of the right ascension- declination type includes the following 
parameters for each recorded spacecraft observation: 

t' , the date and time of observation. The remarks in the main body of this report about re- 
moving reference to the calendar in transforming t' to the relative time, t , also apply here. 

k, the code number for the tracking station reporting the observation 

a^y the observed right ascension, measured in radians eastward from the vernal equinox 

(0 < < 2n). 

the observed declination in radians, measured as positive north of the Equator and as 
negative south of the Equator {-rr/2 < < + rr/2) 

and Wg, the weighting factors corresponding to observations and respectively. This 
information is optional; if not provided, then it is assumed that and wg are each unity. 

The coordinate system employed for the observation data is centered at the tracking station 
on the earth’s surface, and, unlike the system used for recording direction- cosine data, its three 
coordinate axes are parallel to the respective axes of the inertial system. That is, the Z-axis is 
parallel to the earth’s polar axis and the X-Y plane is parallel to the equatorial plane of the earth, 
with the X-axis extending toward the vernal equinox. The Y-axis extends orthogonally to the east 
to form a right-handed system. 

The differential correction process requires the same data to be available as listed in 
the main body of this report, viz., the earth’s flattening coefficient,!, the earth’s rotational 
rate,o), the geodetic longitudes, , of the stations, the geodetic latitudes, 0^, of the stations, 
the altitudes, H , of the stations, the angular distances, , from the vernal equinox to the 
Greenwich meridian at midnight Greenwich mean time for each day in the observational arc, and 
the reference time, 

Computations follow the same scheme given in the main body of this report for the following 
parameters: the auxiliary functions, c and S, the geocentric latitude, the geocentric distance, 
p, of the station, the angular distance, h , between the vernal equinox and the observation meridian 
plane, and the inertial geocentric coordinates Xt, of the station. However, the angle, , 

between the vernal equinox and the tracking station’s X- coordinate axis (measured in the 
observation latitude plane) is zero, since the topocentric and inertial coordinate systems 
are parallel. No rotations are necessary to bring the two systems into coincidence; a single 
translation will suffice. Hence, the relations for the topocentric or local coordinates of the 
satellite are simply: 
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Ym = Y - , 

Zjj = Z - Zj j 

where X, Y, Z are the inertial geocentric coordinates of the satellite predicted by the orbit 
generator. The above simplified relations are obtained from those of the direction-cosine- 
data case by the artificial device of setting = 0 and 6^ - tt/2 in the corresponding equa- 
tions for Y„, Z„ given in the main body of this report (refer to the note at the end of this 
appendix.) 

The computed values of the right ascension and the declination may now be found in terms 
of the local coordinates: 


a 


c 



= arc tan 




1/2 


It is important that the angles and be placed in the proper quadrant for comparison 
purposes with the angles and 8^. In the case of the right ascension, this is done by examin- 
ing the signs of and Y„ separately. The following list presents all possible combinations 


(note that the range for is 0 < 

< 2tt). 




Xm 

> 0, 

Ym 

> 0 : 

0 

7T 

< “= < -J ’ 

X« 

> 0, 

Ym 

< 0 : 

2 

< a < 2 tt 7 

c 

x« 

> 0, 

Ym 

= 0 : 


O 

II 

u 

Xm 

< 0, 

Ym 

> 0 : 

rr 

2 

A 

P 

o 

A 

Xm 

< 0, 

Ym 

< 0 : 

rr 

. _3^ 

< 2 ’ 

Xm 

< 0, 

Ym 

= 0 : 


CL^ ~ IT , 

Xm 

= 0, 

Ym 

> 0 : 


IT 

a = — , 

2 

Xm 

= 0, 

Ym 

< 0 : 


3n 

a = , 

^ 2 

Xm 

= 0, 

Ym 

= 0 : 

a 

c 

indeterminate 


In the case of the declination, the signs of the numerator and denominator of the arctangent 
argument are examined separately. The following list presents all possible combinations (note 
that the range for 8^ is - tt/2 < 8^ < + tt/2). 
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> 0, Z„> 0: 0< S^<f , 


+ >0. Z„ < 0 : --< 8^ < 0 . 


(X5 + y2)‘/^ > 0, Z„ = 0 : 8^ = 0 , 

(Xj + = 0. Z„ > 0 : Sc = f - 

(XS + Yj)*/* = 0, Z„ < 0 : S^ =-| . 

Of course, the case (xj + = Z„ = 0 is not physically possible. 

The observational residuals are now found. 

Aa = ajj - , 


AS = S„ - S . 

0 c 


Here too, care must be exercised. There is one instance where simple subtraction in finding 
the observational residual will yield a misleading result. If one of the right ascensions 
(either observed or computed) is in the first quadrant and very near zero and the other right 
ascension is in the fourth quadrant and very nearly 2 tt , then direct subtraction will provide an 
erroneous result near to 2rr, whereas the intended difference is near to zero. This situation 
can be rectified by the following logical steps: 

If 1^0 ’ ^cl - then Aa = (as above). 

If I - a^l > 7T, then Aa = Sgn (a^ - aj [ 2 tt - \ - a J ] . 

Equivalently, whenever la^ - a^| > 77 , use the following: 

(1) If ^0 ^ then Aa = 2 t 7 - > 0 . 

(2) If ^0 ^ - 2 tt < 0. 

The statistical analysis of the observational residuals follows the procedure given in the 
main body of this report in THE STANDARD DEVIATION OF FIT except that the observational 
residuals are given by Aa. and AS., rather than by AL. and AM.. Hence, the averse residual 
is given by 


n 



The standard deviation of the residuals from their mean value is found from 


cr 



'P [(Aa. - R)^ + (AS. - R)^] 
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The standard deviation of fit is given by 


^6 E ^ 


Modifications are now presented to supplant the material from ANALYTICAL PROCEDURE 
OF DIFFERENTIAL CORRECTION in the main body of this report. 

The first-order Taylor series expansion of the equations of condition may be written 

S T t ^CL 
i=l ^ 


AS 


1 3 S 

i = l ' 


where q. (i = 1, 2, • •• , 6) are the mean or Izsak orbital elements. Expanding the above 
partial derivatives by the chain rule yields 

3^ _ ^ ^ ^ ^ 

3q. 3q. BY^, 3q. 3Z„ 3q. 

^ 

3q. 3X„ 3q. 3Y„ 3q. 3z„ 3qj 

From the equations for and in terms of the local coordinates, we find 


^ - Y„ (X^ + Y^)-^ 


3a 

c 


+ X„ (X^ + Yl) 


3a 

c 


= 0 


3S 

r 


- = - Xm Zm + Z^)-' 
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B 8 

^ = - Y„ Z„ (X2 + yS + zS)-‘ , 

B8 

^ = MXJ + Yi)‘^^ (XS + YS + Z^)-‘ . 


Since the station coordinates X^, Y^, are independent of orbital parameters (and merely 
geodesic functions), the following simple relations hold: 


^ 

^ BZ 

^qi 


The method for calculating the partial derivatives Bx/Bq. , BY/Bq., and B^Bq. is identical to 
that presented in the differential correction scheme in the main body of this report. Then the 
equations of condition are formulated in a precisely analogous manner to that given for the 
direction- cosine data (see THE EQUATIONS OF CONDITION), and there is little need to 
repeat the explicit form of these equations. 

NOTE : The fact that the observational coordinate system is independent of the latitude 

and longitude of the tracking station for right ascension-declination data (as is not the case for 
direction-cosine data) leads to certain possible simplifications in the determination of the 
computed coordinates of the satellite, and 5^. First, recall that the equations for the 
Cartesian inertial coordinates of the observation point are given by 

Xy = yS COS (9 q cos S , 

Y^ = p cos sin 5 , 

= P sin , 

where 8 = (^o)d Here the terms (>^o)d <ij(AT) depend upon the time of the 

observation only, while the term is a function of the location of the observation point. Let 
us denote 


S' = (^o)d + . 


Then we can expand the above equations as 

Xj = p cos Oq cos (S' + k^) 

= p cos 0Q cos cos 8' - cos 0^ sin k^ sin 8' , 
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Yj = p cos 6g sin (8' + \^) 

= p cos 6f, cos sin S' + p cos 0^ sin \g cos 8' , 

Zj - p sin 0g . 


Now denote 


so that 


Xq = p cos 0Q cos , 

Yq = p cos 0g sin \g , 
Zg = p sin 0^ , 


V 


cos 


- sin 

S' 

0 

Yt 

= 

sin 

S' 

cos 

S' 

0 

A 


0 


0 


1 


Xo 

Yo 

Zo 


This represents, in matrix form, the fact that the coordinates Xj, Y^, Zy are obtained from 
Xj, Yq, Zq by a simple rotation about the inertial Z-axis through an angle 8' . Here 8' is the 
angle between the vernal equinox and the Greenwich meridian at observation time. The rec- 
tangular coordinates x„, Y„, z^, obtained directly from the spherical geocentric coordinates 
P> ^o» of tf*® station, represent the Cartesian inertial geocentric coordinates of the track- 
ing station at a time when the Greenwich meridian and the first point of Aries (the vernal 
equinox) coincide. If the coordinates X^, Yg, Zg (dependent upon the station location only) are 
provided as input parameters rather than 0^, and h, then the computations leading up to 
Xj, Y^, Zj are simplified considerably. We need not first compute C, S, 0^, p, and 8. In- 
stead, find s' from parameters relating to the time of observation, and then compute directly 

Xj. = Xg COS S' - Yg sin 8' , 

Yt = Xg sin 8' + Yg cos 8' , 

Zj = Zg . 

Note that this simplified procedure cannot be adopted efficiently with direction-cosine data 
because the rotation matrix involved in computing x^, Y.^, Zj is a function of the geodetic 
latitude, and i/>^, an angular parameter dependent upon the geodetic longitude. 
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